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EXECUTIVE  SUMMARY 

No  accepted  methodologies  exist  to  study  day  to  day  changes  in  the 
characteristics  of  the  circadian  rhythm.  Testing  different  methodologies  using 
data  processing  and  curve-fitting  functions  within  Microsoft  Excel  and  Jandel 
Sigmaplot  is  prohibitively  expensive  in  terms  of  both  time  and  personnel.  Use  of 
the  CIRCAD  program,  described  in  this  report,  dramatically  reduces  the  amount 
of  time  required  for  circadian  data  analyses  and  provides  the  capability  to  quickly 
implement  and  test  new  analytical  methods.  This  Technical  Note  includes  a  brief 
description  of  the  software,  complete  program  listings  in  the  appendices,  and 
sample  input  and  output  including  both  Excel  worksheet  samplers  and  graphs. 
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INTRODUCTION 


PURPOSE 

The  purpose  of  CIRCAD  is  to  automate  the  process  of  analyzing  circadian 
core  temperature  data,  while  providing  the  capability  to  readily  implement  and 
test  new  analytical  methods.  The  analytical  methods  demonstrated  in  the  current 
version  of  CIRCAD  are  variations  of  the  cosinor-fit  method  and  are  used  to 
explore  day  to  day  changes  in  circadian  rhythm  characteristics. 

PROBLEM  STATEMENT 

Various  methodologies  including  spectral  analysis  (e.g.,  Fast  Fourier 
Transform),  autoregressive  modeling  (e.g.,  ARMA  and  MESA),  and  curve-fitting 
approaches  such  as  cosinor  and  Lomb  are  available  for  obtaining  the  circadian 
period  using  core  temperature  data  windows  of  several  days  in  width,  These 
methods  all  assume  “stationary”  data.  For  circadian  data,  stationary  implies  that 
the  period,  phase,  amplitude,  mesor,  zenith,  nadir,  and  acrophase  (Figure  1)  are 
constant  from  one  day  to  the  next.  In  reality,  circadian  data  are  clearly  not 
stationary  (Figure  2).  In  fact,  it  is  the  goal  of  this  project  to  develop  software  to 
aid  in  the  study  of  this  day  to  day  variability. 

Prior  to  the  development  of  the  CIRCAD  software,  implementation  and 
testing  of  analytical  methods  to  obtain  circadian  parameters  on  a  day  to  day 
basis  were  severely  impeded  by  the  time  and  resources  required.  All  data 
manipulation  and  curve  fitting  were  performed  using  Microsoft  Excel  and  Jandel 
Sigmaplot.  Core  temperature  and  actigraphy  data,  including  time  and  date 
marks,  were  copied  to  an  Excel  spreadsheet.  After  copying  the  data  to  Excel, 
technicians  plotted  the  data  to  look  for  outliers  and  valid  data  start  and  stop 
times.  Outliers  and  any  data  collected  outside  the  data  start  and  stop  times  were 
removed  manually.  Core  temperature  and  actigraphy  data  were  then  meshed  to 
a  single  date/time  reference.  Sleep/wake  data  were  entered  manually,  using  the 
same  time/date  reference  as  the  core  temperature  and  actigraphy  data.  Values 
of  0  were  entered  for  all  time  points  in  which  the  subjects  were  awake,  and 
values  of  38.2  were  entered  for  all  time  points  in  which  the  subjects  were  asleep. 
Data  were  then  copied  to  Sigmaplot.  Cosine  curves,  with  a  set  period  length  of 
24  hours,  were  fit  to  each  24-hour  day  (midnight  to  midnight)  within  the  data  set. 
Plots  were  then  generated  showing  the  predicted  versus  the  raw  data.  Scatter 
plots  of  mesor,  period,  amplitude  and  acrophase  versus  day,  relative  to  an  event 
marker,  were  also  generated.  Manipulations  and  curve-fitting,  using  Excel  and 
SigmaPlot,  were  tediously  slow  due  to  the  large  number  of  data  records 
(frequently  exceeding  1 0,000).  Even  after  technicians  became  familiar  with  the 
routine,  analysis  of  each  data  set  averaged  several  man-days.  With  several 
subjects  and  several  data  sets  per  subject,  as  well  as  a  desire  to  investigate 
other  analytical  approaches  besides  cosinor  with  constant  period  length,  it  was 
clear  that  a  more  expedient  process  for  preparing  and  analyzing  data  was 
required. 
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OVERVIEW 

The  CIRCAD  software  first  and  foremost  eliminates  and  speeds  up  much 
of  the  tedium  associated  with  circadian  data  analysis.  Core  temperature  and 
actigraphy  data  must  be  copied  to  the  Excel  file  as  before.  The  technician  must 
enter  into  the  Excel  file  the  time/date  of  all  sleep/wake  events,  the  date  of  the 
reference  event  marker,  and  the  start  and  stop  times  for  valid  data.  The 
technician  then  loads  Matlab  and  runs  the  CIRCAD  program.  It  takes  several 
minutes  for  CIRCAD  to  perform  the  data  manipulation  tasks  of  removing  outliers 
and  meshing  the  core  temperature,  actigraphy,  and  sleep/wake  data  into  a  single 
date/time  reference  that  were  previously  assigned  to  the  technician. 

Once  this  data  manipulation  has  been  accomplished,  CIRCAD  analyzes 
the  data  using  48  different  variations  on  the  cosinor  method.  Other  methods  will 
be  added  in  future  versions  of  the  software.  The  three  types  of  variations  of  the 
cosinor  method  include  (1)  cosine  fits  to  filtered  rather  than  raw  core  temperature 
data,  (2)  alternate  methodologies  for  obtaining  circadian  period,  and  (3) 
alterations  in  the  width  and  starting  point  (in  time)  of  the  data  window  used  in 
performing  the  cosine  curve  fitting.  All  results  are  written  to  the  Excel  file  in  a 
way  that  facilitates  graphing. 


METHODS 

EXPERIMENTAL  DATA  AND  PREPARATION  OF  THE  EXCEL  DATA  FILE 

Three  types  of  data  were  collected.  Core  temperature  was  obtained  using 
a  core  temperature  telemetry  pill  (Human  Technologies  Inc.,  St.  Petersburg,  FL) 
and  a  BCTM2  or  BCTM3  receiver/recorder  (Personal  Electronic  Devices  Inc., 
Wellesley,  MA).  These  data  were  recorded  approximately  once  per  minute. 

Data  were  also  recorded  using  an  actigraphy/skin  temperature  monitor  (Mini- 
Mitter  Company,  Sunriver,  OR).  Actigraphy  data  consist  of  the  counts  of  body 
movements  in  a  1 -minute  period.  In  a  few  cases,  skin  temperature  data  at  2-3 
locations,  were  also  recorded.  Finally,  sleep/wake  data  were  recorded  manually, 
on  a  log  sheet,  kept  by  the  test  subjects.  Each  type  of  data  was  associated  with 
its  own  time  and  date  marks.  Data  were  typically  recorded  for  a  2-10  day  period 
from  a  single  subject  for  each  data  set. 

Core  temperature  and  actigraphy  data  were  first  downloaded  from  the 
recording  device  to  an  Excel  worksheet  on  a  personal  computer.  Separate 
date/time  references  were  maintained  for  each  recording  device  (i.e.,  BCTM  and 
MiniMitter).  Sleep/wake  data  were  manually  entered  into  the  worksheet.  Sleep 
or  wake  onset  date  and  time  were  entered  as  separate  columns;  a  thirdO  column 
was  used  to  contain  the  sleep/wake  code  (“S”  for  sleep  onset,  or  “W”  for 
awakening).  The  technician  followed  strict  rules  in  naming  variables  within  the 
worksheet  and  ensured  that  there  were  no  blank  values  in  any  of  the  time  or  date 
columns.  A  file  naming  system  was  developed  that  includes  codes  to  identify  the 
subject,  the  data  set,  and  the  Julian  date  and  year  at  the  start  of  data  collection. 
The  worksheet  containing  the  raw  BCTM  and  MiniMitter  data  was  given  the  same 
name.  This  was  done  because  the  functions  for  dynamic  data  exchange  (DDE) 
in  Matlab  make  connections  to  worksheets  of  a  specified  name;  these  DDE 
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functions  ignore  file  name  specifications.  Thus,  if  the  same  worksheet  name 
were  to  be  used  in  2  different  files,  Matlab  would  read  data  from  the  worksheet  of 
whichever  file  was  active  at  the  time.  The  only  way  to  ensure  that  the  correct 
worksheet  is  read  by  CIRCAD  is  to  use  worksheet  names  that  are  identical  to  the 
file  names. 

Once  the  raw  data  worksheet  was  generated  and  properly  named,  a  new 
worksheet  named  Genlnfo  was  created.  The  reference  event  date  and  the  data 
start  and  stop  times  are  then  entered  into  the  Genlnfo  worksheet.  Two  other 
worksheets,  MATLAB  Summary  Results  and  MATLAB  Results  for  Graphs, 
were  also  created  to  which  results  can  be  written.  A  sample  portion  of  an  Excel 
worksheet  is  shown  in  Figure  3. 

DEVELOPMENT/DEPLOYMENT  ENVIRONMENT 

All  programs  are  implemented  in  Matlab  script  language  (Mathworks,  Inc., 
Natick,  MA)  and  run  within  the  Matlab  environment.  For  the  purposes  of 
modifying  and  enhancing  the  software  in  the  future,  functions  were  encapsulated 
to  the  greatest  practical  extent.  In  Matlab,  external  functions  are  required  to  be 
saved  in  a  Matlab  script  file  (m-file)  of  the  same  name.  A  list  of  all  the  m-files  that 
make  up  the  CIRCAD  software  is  provided  in  Appendix  A. 


THE  MAIN  FUNCTION  (CIRCAD) 

The  “main”  function,  used  to  call  all  other  functions,  is  called  circad.  circad  is 
located  within  the  Matlab  m-file  of  the  same  name  (see  Appendix  B).  circad  calls 
the  function  ReadCircadData,  which  reads  the  core  temperature,  actigraphy,  and 
sleep/wake  data  from  the  Excel  file.  Meshing  the  data  into  a  single  date/time 
reference  with  a  constant  sampling  interval  is  accomplished  using  the  function 
MeshData.  Sleep/Wake  data  are  transformed  to  a  time-based  vector  using  the 
same  date/time  reference  as  the  core  temperature  and  actigraphy  data  in 
TransformSIeepWake.  “S”  values  are  recoded  as  38.2,  and  “W”  values  as  0  in 
TransformSIeepWakeData  as  well.  Removal  of  outliers  is  accomplished  in 
RemoveOutliers.  Filtering  of  data  is  performed  in  the  function  FFT_Filtering, 
which  uses  Matlab’s  built-in  fast  Fourier  transform  function.  The  function 
RunAnalyses  sets  parameter  values  to  be  used  in  the  model  fitting  procedures, 
calls  the  DoCosinor function,  and  prints  and  saves  the  results.  DoCosinor is 
responsible  for  fitting  a  cosine  curve  to  specified  data  intervals.  All  of  these 
functions  are  listed  in  the  appendices  along  with  any  supporting  functions  and 
are  described  in  more  detail  below. 


READING  CIRCADIAN  DATA  ( READCIRCADDATA ) 

ReadCircadData  (Appendix  C)  is  responsible  for  reading  data  from  the 
Excel  file.  This  is  accomplished  using  Matlab’s  built-in  functions  for  dynamic  data 
exchange  (DDE).  At  the  start  of  the  project,  data  existed  in  Excel  worksheets  in 
many  different  formats.  To  minimize  the  amount  of  reformatting  required  and  to 
take  advantage  of  data  manipulation  that  had  already  been  performed  manually, 
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ReadCircadData  was  written  to  provide  a  certain  amount  of  flexibility.  Although  it 
was  necessary  to  insist  on  a  standard  set  of  variable  names,  it  was  not 
necessary  to  insist  on  standard  locations  for  each  variable  within  the  data  file.  It 
was  also  possible  to  accommodate  the  three  types  of  “date”  formats  that  had 
been  used;  “date”  could  be  entered  either  as  a  day  of  the  month  (Day),  a  date 
string  (Date)  or  the  number  of  elapsed  days  from  the  start  of  data  collection 
(EIDays).  Recognized  variable  names  include  the  following:  BCTM-Time,  BCTM- 
Date,  BCTM-Day,  BCTM-EIDays,  Pill  Temp,  Mini-Time,  Mini-Date,  Mini-Day, 
Mini-EIDays,  Twrist,  Tother,  Activity,  Sleep-Time,  Sleep-Date,  Sleep-EIDays,  and 
Sleep/Wake.  Variable  names  may  occur  anywhere  within  the  first  30  columns 
and  15  rows  of  the  worksheet.  It  is  assumed  that  whatever  falls  below  the 
variable  name  (in  the  same  column)  is  the  data  for  that  variable.  The  only 
mandatory  columns  are  BCTM-Time,  and  one  of  the  BCTM  "date"  options  (i.e., 
BCTM-Date,  BCTM-Day,  or  BCTM-EIDays).  If  actigraphy  or  sleep/wake  data 
occur  without  corresponding  date  and  time  columns,  they  are  assumed  to  occur 
at  the  same  time/date  points  as  the  BCTM  data  on  the  same  row.  ~~ 

ReadCircadData  initially  reads  data  from  an  Excel  file.  However,  because 
reading  from  Excel  files  is  a  slow  process,  the  data  are  saved  to  a  Matlab  binary 
file  (mat-file)  for  subsequent  sessions  (mat-files  can  be  read  much  faster  than  the 
corresponding  Excel  files.  ReadCircadData  will  read  the  mat  file  if  it  exists. 
Unless  otherwise  specified,  all  of  the  discussions  that  follow  pertain  to  reading 
data  from  Excel  files. 

ReadCircadData  performs  the  following  tasks:  (1)  get  the  Excel  file  name 
using  GetFileName,  (2)  check  to  see  if  a  mat-file  (Matlab  binary  file)  exists  using 
the  function  DecideBetweenExcelAndMatFile,  (3)  read  the  Excel  or  mat-file  data 
(ReadExcelData  or  ReadMatData),  and  (4)  if  the  data  came  from  an  Excel  file, 
save  the  raw  data  as  a  mat  file  using  the  function  Save2MatFile. 

At  the  conclusion  of  ReadCircadData,  a  structure  X  exists  that  contains  all 
the  data  read  from  the  Excel  or  mat  file.  X  currently  includes  the  following  fields: 


PathName 

Data  path  name  (string) 

FileName 

■  III  1 1 1  1 II  III  1 1— 

FileHandle 

File  name  without  the  extension  (string) 

BCTM_DateNum 

Time  and  Date  from  BCTM  as  a  date  number 
(numeric  array) 

BCTMJHour 

Time  and  Date  from  BCTM  as  the  elapsed 
hours  since  midnight  of  the  first  day  (numeric 
array) 

PillTemp 

Pill  temperature  data  (numeric  array) 

Mini_DateNum 

MiniJHour 

Time  and  Date  from  MiniMitter  as  the 
elapsed  hours  since  midnight  of  the  first  day 
(numeric  array) 

Twrist 

Wrist  temperature  data  (numeric  array) 

Tother 

Auxiliary  skin  temperature  data  from  Mini- 
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Mitter  (numeric  array) 

Activity 

Actigraphy  data  (numeric  array) 

Sleep_DateNum 

Time  and  Date  from  the  sleep/wake  log  as  a 
date  number  (numeric  array) 

Sleep_Hour 

Time  and  Date  from  the  sleep/wake  log  as 
the  elapsed  hours  since  midnight  of  the  first 
day  (numeric  array) 

SleepWake 

Sleep/Wake  onset  from  sleep/wake  log 
(numeric  or  string  array) 

StartDateNum 

Start  Date  for  data  collection  as  a  date 
number 

EventDateNum 

Date  of  the  reference  event  converted  to  a 
date  number 

DataStartDateNum 

Start  of  “good”  data  as  a  date  number 

DataStopDateNum 

End  of  “good”  data  as  a  date  number 

MESHING  DATA  INTO  A  SINGLE  TIME/DATE  REFERENCE  (MESHDATA) 

The  process  of  meshing  the  BCTM,  MiniMitter  and  Sleep/Wake  data  into  a 
single  Date/Time  reference  is  straightforward,  but  requires  a  large  amount  of  time 
if  performed  manually.  The  function  MeshData  (Appendix  D)  is  found  in  a  Matlab 
m-file  of  the  same  name.  The  first  step  in  MeshData  is  to  copy  all  of  the  non¬ 
time-based  fields  in  X  (i.e.,  PathName,  FileName,  FileHandle,  StartDateNum, 
EventDateNum,  DataStartDateNum,  and  DataStopDateNum)  to  a  new  structure, 
Y.  This  is  accomplished  in  the  function  CopyNonMeshedFields2Y.  The  next 
step  is  to  compute  the  appropriate  start  and  stop  times  for  valid  data  using  the 
function  GetDataStartStopTimes.  The  start  time  is  either  the  first  valid  date/time 
in  any  of  the  device  fields,  or  the  data  start  time  as  entered  into  the  Genlnfo 
worksheet,  whichever  is  later.  The  stop  time  is  either  the  last  valid  date/time  in 
any  of  the  device  fields,  or  the  data  stop  time  as  entered  into  the  Genlnfo 
worksheet,  whichever  is  earlier.  Next,  a  new  array  field  is  created  in  structure  Y, 
named  Y.Time,  to  hold  time  (in  hours).  The  first  value  in  Y.Time  is  the  start  time, 
the  last  value  is  the  stop  time,  and  the  interval  between  data  points  (Y.dTime)  is 
one  (1)  minute  or  1/60  hours.  All  of  the  time/date  fields  are  rounded  to  the 
nearest  minute.  Finally,  all  of  the  time-based  arrays  are  mapped  into  the  new 
time  reference  using  the  function  DoTheMapping. 

The  function  DoTheMapping  and  the  function  it  calls,  Map,  both  take 
advantage  of  a  special  Matlab  feature  that  allows  array  assignments  to  be  made 
using  arrays  of  index  values.  For  example,  if  x=[1 1  22  33  44  55]  and  i=[2  5],  then 
the  statement  y=x(i)  will  result  in  y=[22  55].  It  is  also  possible  to  use  the 
statement  y(i)=x(4:5),  which  would  result  in  y=[NaN  44  NaN  NaN  55],  where  NaN 
is  an  empty  placeholder  in  Matlab.  NaN  stands  for  "not  a  number". 

For  any  of  the  devices  (i.e.,  BCTM,  MiniMitter  or  Sleep/Wake  log),  the 
difference  in  minutes  between  the  device_Hour  value  and  the  start  time  will 
produce  the  position  of  the  device's  data  values  with  respect  to  the  new  time 
reference.  Adding  1  to  these  values  gives  the  appropriate  array  index.  These 
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array  indices  are  named  BCTM.i,  Mini.i,  and  Sleep.i.  Do TheMapping  calls  the 
function  Map  tor  each  of  the  time-based  data  arrays  (i.e.,  PillTemp,  Twrist, 
Tother,  Activity,  and  SleepWake).  When  Map  is  called,  the  data  array  is  passed 
into  the  variable  Data,  and  the  array  indices  are  passed  into  the  variable  iData. 

N  is  the  size  of  the  Y.Time  array,  and  NumericYN  tells  Map  whether  the  data 
array  contains  numeric  data  (1  for  numeric,  0  for  string). 

The  first  thing  Map  does  is  check  to  see  if  Data  is  empty.  If  it  is,  then  the 
mapped  variable  must  also  be  empty.  If  Data  is  not  empty,  then  the  mapping 
process  is  started.  First,  Data  and  iData  are  trimmed  to  remove  elements 
corresponding  to  data  collected  before  the  start  time  or  after  the  stop  time.  This 
is  done  by  identifying  values  in  iData  which  are  greater  than  or  equal  to  1  and 
less  than  or  equal  to  N.  This  valid  interval  for  Data,  called  iGoodData,  is  used  to 
trim  the  arrays:  Data=Data(iGoodData)  and  iData=ix(iGoodData).  The  array  that 
Data  is  to  be  mapped  into  (DataMapped)  is  first  filled  with  empty  placeholders 
(numeric  arrays  only).  The  empty  placeholder  for  numeric  arrays  in  Matlab  is  the 
value  NaN.  The  mapping  process  itself  is  accomplished  using  the  command 
DataMapped=Data(iData).  When  Matlab  assigns  array  values  using  a  variable 
index,  the  size  of  the  new  array  is  set  to  the  largest  value  in  the  array  index.  In 
this  case,  the  size  of  DataMapped  is  set  to  the  largest  value  in  iData.  The  largest 
value  in  iData  is  not  necessarily  equal  to  N,  resulting  in  DataMapped  being 
smaller  than  the  time  reference  array  (Y.Time).  To  get  around  this,  a  series  of 
statements  in  Map  beginning  with  "if  length(DataMapped)<N"  assigns  a  value  to 
the  Nth  value  of  DataMapped  (NaN  for  numeric  arrays  and  the  last  "good"  string 
value  for  string  arrays)  which  forces  the  size  of  the  mapped  array  to  be  equal  to 
N.  The  final  meshed  data  are  named  Y.Time,  Y.Traw,  Y.Twrist,  Y.Tother, 
Y.Activity,  and  Y.SIeepWake.  All  of  these  arrays  except  Y.SIeepWake  are 
written  to  the  Excel  worksheet  named  MATLAB  Results  for  Graphs.  Upon 
returning  to  the  calling  function,  circad,  X  is  reassigned  the  values  in  Y. 


TRANSFORM  SLEEP/WAKE  DATA  ( TRA NSFORMSL E EP WA KE) 

The  function  TransformSIeepWake  (Appendix  E)  takes  sleep/wake  event 
data  (X.SIeepWake),  where  the  event  "sleep"  is  coded  either  with  the  numerical 
value  of  38.2  or  the  letter  “S”,  and  “wake"  is  coded  either  with  the  numerical 
value  of  0  or  the  letter  W.  For  string  data  (i.e.,  letter  codes), 
TransformSIeepWake  converts  these  data  to  a  numerical  array  in  which  S  and  W 
are  represented  by  their  corresponding  ASCII  codes  (83  and  87,  respectively) 
and  missing  values  are  coded  zero  (0).  If  the  first  value  in  the  array  is  missing 
(likely  because  the  data  have  already  been  mapped  to  the  common  date/time 
reference),  it  is  coded  with  the  value  87  (assuming  that  the  subject  is  awake  at 
the  start  of  data  collection).  Other  missing  values  are  filled  in  by  dragging  the 
previous  value  along.  After  all  missing  values  have  been  filled  in,  all  values  of  83 
are  converted  to  38.2  (a  convenient  value  for  representing  sleep  periods  on  core 
temperature  plots),  and  values  of  87  are  converted  to  the  value  0.  Transformed 
SleepWake  data  are  written  to  the  Excel  worksheet  named  MATLAB  Results  for 
Graphs  prior  to  returning  to  the  calling  function. 
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REMOVING  OUTLIERS  (REMOVEOUTLIERS) 

The  function  RemoveOutliers  (Appendix  E)  takes  the  array  X.Traw  and 
copies  it  to  X.Tclean.  Core  temperature  values  in  X.Tclean  that  are  strings,  or 
that  are  less  than  35  or  greater  than  40  are  set  equal  to  NaN.  Cleaned-up  core 
temperature  data  (X.Tclean)  are  written  to  the  Excel  worksheet  named  MATLAB 
Results  for  Graphs  prior  to  returning  to  the  calling  function. 

FILTERING  USING  FFT  (FFT_FILTERING) 

Filtering  is  performed  using  the  Fast  Fourier  Transform  (FFT).  The 
function  FFT_Filtering  (Appendix  F)  uses  Matlab’s  built-in  FFT  ( fft )  and  Inverse 
FFT  ( ifft )  functions.  FFT_Filtering  starts  with  the  cleaned-up  core  temperature 
vector,  X.Tclean.  To  use  fft  and  ifft,  it  is  necessary  for  data  to  be  equally  spaced 
in  time.  Because  the  time  array  (X.Time)  is  already  equally-spaced  at  1  -minute 
intervals,  this  requires  that  all  missing  temperature  values  (i.e.,  values  of  NaN  in 
X.Tclean)  be  assigned  some  kind  of  estimated  temperature  value. 

Filling  in  missing  temperature  values  is  accomplished  in  the  function 
FillTemperatureBlanks.  For  FFT  analysis,  missing  value  estimation  may  consist 
of  either  setting  the  missing  data  points  to  the  average  value  of  the  entire  data 
series,  or  of  setting  any  missing  data  point  equal  to  the  previous,  good  data 
point1.  The  latter  approach  is  used  in  FillTemperatureBlanks.  It  is  possible  that 
the  first  data  point(s)  in  the  data  series  may  also  contain  missing  data.  In  this 
case,  starting  with  the  first  “good”  data  point  and  moving  backwards  in  the  series, 
any  missing  data  point  is  set  equal  to  the  next  good  data  point  in  the  series. 

Once  a  valid  core  temperature  vector  is  obtained  (named  X.T),  the  FFT  is 
performed,  producing  the  vector  FFTofT.  FFT  frequency  (freq),  period,  and 
magnitude  vectors  are  also  computed. 

Filtering  is  accomplished  by  selecting  a  frequency  cut-off  point,  and  setting 
all  elements  of  FFTofT  that  correspond  to  frequencies  greater  than  the  cut-off 
value  to  zero.  The  cut-off  frequencies  used  in  this  version  of  CIRCAD  are  0.25, 
0.50,  and  1 .00  hours'1,  corresponding  to  cut-off  periods  of  4,  2,  and  1  hours, 
respectively.  Taking  the  inverse  of  the  filtered  FFTofT  vector  produces  the 
“filtered”  vectors  of  core  temperature  with  respect  to  time  (i.e.,  X.Tfiltl ,  X.Tfilt2, 
X.Tfilt4,  which  correspond  to  cut-off  periods  of  1, 2,  and  4  hours,  respectively). 

The  final  step  in  the  FFT_Filtering  function  is  to  plot  the  results  and  send 
the  filtered  arrays  (X.Tfiltl ,  X.Tfilt2,  and  X.Tfilt3)  to  the  Excel  worksheet  MATLAB 
Results  for  Graphs. 

RUN  ANALYSES 

All  analyses  are  performed  by  calling  the  RunAnalyses  function  (Appendix 
B)  from  within  circad.  RunAnalyses  first  selects  one  of  the  48  different  conditions 
to  run.  Each  condition  includes: 

1 .  The  variable  (Var)  to  fit  (i.e.,  X.T,  X.Tfiltl ,  X.Tfilt2,  or  X.Tfilt4); 

2.  The  method  (PeriodStr)  for  obtaining  the  circadian  period  (i.e,  “free”  in 
which  the  cosine-fit  procedure  estimates  the  period,  “24-hour”,  in  which 
a  constant  period  of  24  hours  is  used,  or  “Lomb”  in  which  the  function 
LombPeriod  (Appendix  G)  is  used  to  compute  the  circadian  period); 
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3.  The  data  interval  (Interval)  to  use  (i.e.,  whole  curve,  individual  24-hour 
(midnight  to  midnight)  days,  individual  24-hour  (noon  to  noon)  days,  or 
36-hour  (6pm  to  6am)  segments). 

RunAnalyses  then  calls  the  appropriate  analytical  function,  passing  the 
appropriate  time  and  temperature  vectors,  as  well  as  a  numeric  code 
representing  the  period,  in  the  function  call.  In  the  current  version  of  CIRCAD, 
there  is  only  one  analytical  function  named  DoCosinor,  which  performs  a  cosine 
fit  to  the  data. 

The  result  of  calling  the  analytical  function  (e.g.,  DoCosinoi)  is  to  return  a 
data  structure  containing  the  results  to  the  variable  IntervalResults. 
IntervalResults  may  describe  either  a  fit  to  the  entire  data  set,  or  to  a  24-  or  36- 
hour  segment  of  the  data  set.  For  the  case  in  which  the  interval  is  the  entire  data 
set,  IntervalResults  is  copied  directly  into  a  structure  named  Results.  In  the  case 
where  the  intervals  are  24-  or  36-hours,  IntervalResults  structures  are  appended 
together  into  the  Results  structure  to  cover  the  entire  period  of  data  collection. 

For  example,  if  there  are  5  days  of  data,  then  DoCosinor  is  run  once  for  each  of 
the  5  days.  The  resulting  IntervalResults  structures  are  appended  together  to 
form  the  structure  named  Results.  Although  the  Results  structure  contains  the 
same  fields  as  IntervalResults,  the  fields  are  generally  arrays  rather  than  scalars 
(one  array  element  for  each  day).  Also,  the  Tpred  and  t  arrays  would  contain 
data  from  all  5  days.  Note  that  in  the  case  of  36-hour  data  windows,  values  of 
Tpred  and  t  from  midnight  to  midnight  are  appended  together;  the  6-hour  data 
segments  preceding  and  succeeding  the  24-hour  period  are  not  saved. 

Once  the  results  are  appended  together,  predicted  values  are  plotted 
along  with  whatever  data  were  used  to  obtain  the  predicted  values  (i.e.,  “cleaned- 
up”  or  filtered).  Predicted  values  are  also  plotted  against  the  unfiltered  data. 
These  plots  are  automatically  saved  as  .jpg  files.  Predicted  values  are  also 
written  to  the  Excel  worksheet  named  MATLAB  Results  for  Graphs  for  later 
plotting.  Values  of  summary  variables  (e.g.,  mesor,  amplitude,  period,  etc.)  are 
written  to  the  Excel  worksheet  named  MATLAB  Summary  Results. 


COSINOR  ANALYSES  (DOCOSINOR) 

DoCosinor  is  called  from  RunAnalyses  to  perform  a  cosine  fit  to  the  data 
passed  to  it.  DoCosinor  is  listed  in  Appendix  H  and  is  responsible  for  computing 
the  following  fields  of  the  structure  named  cosinor: 


Tpred 

Predicted  temperature  using  cosinor  model  and 
fit  parameters 

t 

The  time  interval  used 

mesor 

The  mesor  (y-offset,  estimated  mean 
temperature)  as  computed  by  the  cosinor 
regression. 

amplitude 

The  amplitude  as  computed  by  the  cosinor 
regression 

period 

The  period  as  computed  by  the  cosinor 
regression  (if  applicable),  24,  or  the  Lomb 
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period. 

phase 

The  phase  as  computed  by  the  cosinor 
regression  (for  example,  -8.75) 

gphase 

The  phase  relative  to  the  sampling  interval  (for 
example,  39.35) 

clock_phase 

The  phase  relative  to  clock  time  (for  example, 
15.25). 

1|H  ■  ||  ||— 

The  phase  as  a  text  string  (for  example,  15:15). 

R 

Correlation  coefficient  (R-value)  for  the 
regression 

N 

Number  of  data  points  used  in  the  regression 

MaxData 

Maximum  of  the  data  values 

TimeOfMaxData 

Time  (clock  hours)  corresponding  to  the 
maximum  of  the  data  values 

MinData 

Minimum  of  the  data  values 

TimeOfMinData 

Time  (clock  hours)  corresponding  to  the 
minimum  of  the  data  values 

MaxPred 

Maximum  of  the  predicted  values 

TimeOfMaxPred 

Time  (clock  hours)  corresponding  to  the 
maximum  of  the  predicted  values 

MinPred 

Minimum  of  the  predicted  values 

TimeOfMinPred 

Time  (clock  hours)  corresponding  to  the 
minimum  of  the  predicted  values 

Mesor,  amplitude,  period,  and  the  various  acrophase  parameters  (phase, 
gphase,  clock_phase,  and  phase_clockstring)  are  obtained  by  fitting  a  cosine 
function  to  the  time  and  Temperature  data  that  are  passed  to  the  DoCosinor 
function.  To  obtain  values  for  these  fields,  DoCosinor  first  assigns  initial  values 
to  the  3  or  4  parameters  (p)  that  will  be  fit.  The  number  of  parameters  depends 
on  whether  the  period  is  to  be  computed  during  the  cosine  fit  procedure 
(Period==0),  or  whether  period  is  set  at  24-hours  (Period==24)  or  computed 
using  the  Lomb  analysis  method.  The  first  three  initial  parameter  values 
correspond  to  the  mesor  (°C),  amplitude  (°C)  and  acrophase  (radians),  and  the 
intial  values  are  p(1  )=40,  p(2)=0.5,  and  p(3)=1 ,  respectively.  If  the  period  is  to  be 
computed  using  the  cosine  fit  procedure,  then  a  fourth  parameter  is  used  which 
corresponds  to  the  frequency  (radians/hour)  with  an  initial  value  of  p(4)=0.25. 
These  initial  values  were  selected  arbitrarily  and  have  been  found  to  produce 
satisfactory  fits  for  all  data  analyzed  to  date. 

Once  initial  values  have  been  specified,  DoCosinor  calls  either  sumsqerror 
or  sumsqerrorFP,  depending  on  whether  the  period  is  to  be  computed  or  has 
already  been  set,  respectively.  The  sumsqerror  and  sumsqerrorFP  functions 
minimize  the  prediction  error  of  the  function  CosinorTpred,  which  uses  the  cosine 
function  to  predict  core  temperature: 

Tpred=p(1)+p(2)*cos(p(4)*t+p(3)) 

where  2*pi/period  is  substituted  for  p(4)  when  the  period  is  not  to  be  computed 
using  the  cosinor  fit  procedure. 

The  values  of  the  cosinor  fields  are  computed  after  the  completion  of  the 
parameter  optimization  step: 
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•  Tpred  is  computed  by  calling  the  CosinorTpred  function  using  the  optimized 
parameter  values 

•  t  is  copied  from  the  time  vector  that  is  passed  to  DoCosinor 

•  N  is  set  to  the  length  of  the  t  vector 

•  R,  the  correlation  coefficient  between  T  and  Tpred,  is  computed  using 
Matlab’s  built-in  function  named  corrcoef 

•  mesor  and  amplitude  correspond  directly  to  the  first  and  second  optimized 
parameter  values,  p(1)  and  p(2),  respectively 

•  period  is  given  by  p(4),  the  frequency  parameter,  by  taking  the  reciprocal  of 
p(4),  then  converting  from  radians  to  cycles:  period  =  2n/p(4). 

•  phase  (“acrophase”,  in  hours)  is  computed  using  p(3),  the  phase  in  radians, 
and  the  period  (in  hours/cycle,  computed  in  the  step  above): 

phase  =  -[p(3)/(2*7c)]*period 

For  ease  in  viewing  the  location  of  the  phase  on  the  composite  graphs 
(with  time  in  hours  on  the  x-axis),  “gphase”  is  computed  by  repeatedly  adding  the 
period  length  to  phase  (gphase  =  phase  +  period  +  period  +...)  until  gphase  is 
greater  than  the  first  time  value  t(1)  in  the  interval.  “clock_phase”  (i.e., 
acrophase  in  hours  after  the  previous  midnight)  and  “phase_clockstring”  (i.e., 
clock_phase  as  a  time  string)  are  computed  by  repeatedly  subtracting  24  hours 
from  phase  (clock_phase  =  phase  -24  -24  -...)  until  a  valid  clock  time  (i.e., 
between  0  and  24  hours)  is  obtained. 

Matlab’s  built  in  max  and  min  functions  are  applied  both  to  the  data  values 
(T)  and  the  predicted  values  (Tpred)  to  obtain  both  the  max  and  min  values,  and 
also  the  array  indices  corresponding  to  the  max  and  min  values.  From  the 
indices,  it  is  possible  to  obtain  the  corresponding  time  of  the  max  and  min  values. 
The  final  time  values  are  computed  by  repeatedly  subtracting  24  hours  from  the 
time  value  (e.g.,  TimeOfMax  =  TimeOfMax  -24  -24  -...)  until  a  valid  clock  time 
(i.e.,  between  0  and  24  hours)  is  obtained. 

MISCELLANEOUS  FUNCTIONS 

Utility  functions  that  are  called  by  more  than  one  function  or  that  are 
utilized  by  other  applications  are  contained  in  their  own  m-files.  These  m-files 
are  listed  in  Appendix  I. 


RESULTS 


SAMPLE  INPUT/OUTPUT 

Sample  output  from  CIRCAD  are  shown  in  Figures  4  through  8.  Three 
types  of  figures  are  produced  in  MATLAB  during  CIRCAD  program  execution. 
The  first  shows  the  results  of  the  FFT_Filtering  process  (Figure  4).  The  topmost 
plot  is  of  the  raw  data.  The  next  plot  shows  the  FFT-generated  periodogram. 

The  final  three  plots  show  the  filtered  data  (frequency  cut  offs  set  to  1 ,  0.5,  and 
0.25  hours'1,  respectively).  The  next  type  of  figure  (Figure  5)  shows  the  results  of 
one  of  the  48  different  variations  of  the  cosinor  method.  In  the  example  shown,  a 
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36-hour  data  window  was  used,  the  cosinor  was  used  to  fit  the  circadian  period, 
and  filtered  data  (f>0.5  h"1)  was  used.  For  those  variations  in  which  the  Lomb 
method  was  used  to  obtain  the  dominant  (i.e.,  circadian)  period,  a  Lomb 
periodogram  is  also  displayed.  The  example  Lomb  periodogram  in  Figure  6 
corresponds  to  a  variation  using  a  36-hour  data  window  using  cleaned-up  (not 
filtered)  data. 

CIRCAD  also  writes  results  back  to  the  Excel  data  file  so  that  publication- 
quality  graphics  can  be  generated  using  SigmaPlot.  The  time-based  data, 
including  the  raw,  cleaned  up,  and  filtered  core  temperature  data  as  well  as  all  of 
the  cosinor  fits  are  written  to  the  worksheet  named  MATLAB  Results  for 
Graphs.  An  example  of  the  part  of  the  worksheet  containing  the  raw,  cleaned 
up,  and  filtered  core  temperature  data  is  shown  in  Figure  7a.  Figure  7b  shows 
the  results  of  one  of  the  48  different  variations  of  the  cosinor  method  in  which  a 
36-hour  data  window  was  used,  the  cosinor  was  used  to  fit  the  circadian  period, 
and  cleaned-up  data  (not  filtered)  was  used  for  the  fit. 

In  addition  to  the  time-based  results,  CIRCAD  also  writes  summary  data, 
containing  the  parameter  values  for  each  of  the  cosinor  fits,  to  the  worksheet 
named  MATLAB  Summary  Results.  An  example  from  a  MATLAB  Summary 
Results  worksheet  showing  the  results  of  a  single  variation  of  the  cosinor  fit 
method  is  provided  in  Figure  8. 


TIME  SAVINGS 

Manually  performing  the  data  manipulation,  curve  fitting  and  plotting  using 
Excel  and  Sigmaplot  required  approximate  3Vz  days  for  one  cosinor  fit  (i.e.,  a 
cosinor  using  a  constant  period  of  24  hours  and  midnight  to  midnight  data 
windows  on  raw  data).  Using  the  CIRCAD  software,  this  process  was  shortened 
considerably.  It  now  requires  approximately  30  minutes  to  read  and  process  the 
data,  and  to  perform  48  variations  of  the  cosinor  method  (i.e.,  variations  in 
circadian  period,  data  windows,  and  data  filtering).  It  requires  an  additional  20 
minutes  to  plot  the  results,  which  is  performed  using  Sigmaplot  templates. 

DISCUSSION:  PROBLEMS/LIMITATIONS  AND  FUTURE 

DIRECTIONS 

The  current  version  of  CIRCAD  utilizes  only  the  cosinor  method  for 
analyzing  core  temperature  data.  We  intend  to  explore  other  signal  processing 
and  curve-fit  methods  for  obtaining  circadian  parameters  in  future  versions  of  the 
software.  Also,  we  would  like  to  apply  these  methods  to  activity  and  skin 
temperature  data  as  well  as  body  core  temperature  data. 

Although  the  CIRCAD  software  significantly  reduces  the  amount  of 
manual  labor  involved  in  analyzing  body  core  temperature,  there  is  still  a  large 
amount  of  tedious  labor  (approximately  30  minutes  per  data  set).  We  would  like 
to  reduce  this  as  far  as  is  practical.  This  may  include  obtaining  sleep/wake  data 
directly  from  the  actigraphy  monitor,  modifying  the  way  data  are  read  and  written 
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to  the  Excel  worksheets,  and  creating  a  method  for  performing  sequential 
analyses  on  multiple  data  sets  (i.e.,  enabling  a  "batch"  mode). 

Certain  limitations  are  also  imposed  by  current  “bugs”  and  limitations  in 
the  DDE  functions  in  the  Matlab  programming  environment.  Currently,  it  is  not 
possible  to  automatically  start  up  Excel  and  open  the  desired  input  data  file;  the 
user  must  perform  these  operations  manually.  It  is  not  possible  to  read  from 
more  than  one  data  file,  or  to  read  from  one  file  and  place  results  in  another 
(although  it  is  possible  to  work  with  multiple  worksheets  within  a  single  Excel  file). 
Because  the  file  is  normally  re-saved  after  results  are  written  to  it  from  CIRCAD, 
the  number  of  instances  in  which  it  is  possible  to  save  time  by  reading  from  a 
mat-file  is  severely  reduced.  Finally,  as  was  mentioned  previously,  the  DDE 
functions  do  not  pay  attention  to  workbook  (i.e.,  file)  names,  so  input  worksheet 
names  must  be  unique  to  prevent  mistakes  in  reading  the  data.  The  next  version 
of  the  CIRCAD  software  will  likely  utilize  ActiveX  (as  opposed  to  DDE)  methods 
to  communicate  with  Excel.  Matlab’s  ActiveX  functions  have  so  far  shown  to 
overcome  many  of  the  limitations  associated  with  DDE  in  terms  of  reading  data 
from  Excel.  Also,  ActiveX  is  a  newer  technology  that  is  rapidly  replacing  DDE  as 
the  standard. 


13 


REFERENCES 


1.  Press,  W.H.,  B.P.  Flannery,  S.A.  Teukolsky,  and  W.T.  Vetterling.  Numerical 
Recipes  in  C:  The  Art  of  Scientific  Computing.  Cambridge  University  Press, 
New  York,  1992. 

2.  Coyne,  M.D.  and  L.A.  Stephenson.  Personal  Communication,  1999. 


14 


Figure  1.  Parameters  of  a  Circadian  (approximately  24-hour)  Rhythm. 


Figure  2.  Plot  demonstrating  day  to  day  variability  in  the  circadian  core 
temperature  rhythm  (from  M.  Coyne  and  L.  Stephenson,  unpublished2). 
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Figure  3.  A  sample  from  an  Excel  data  worksheet. 
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Figure  4.  Example  plot  produced  by  the  FFT_Filtering  function.  The  first 
graph  is  of  core  temperature  data  (X.T)  versus  time.  The  second  graph  is 
the  FFT  periodogram.  The  last  three  graphs  show  the  effects  of  filtering 
the  core  temperature  data  using  cutoff  periods  of  1-hour,  2-hours,  and  4- 
hours,  respectively. 
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Figure  5.  Example  MATLAB-generated  plots  showing  the  cosinor  fit 
against  the  cleaned-up  raw  data  (top)  and  also  against  the  variable  that  was 
fit,  in  this  case  filtered  (cutoff  frequency  =  0.5  hours'1)  core  temperature 
data. 
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Figure  6.  Example  of  a  Lomb  Periodogram  produced  by  CIRCAD. 
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Figure  7a.  Example  output  from  CIRCAD.  This  figure  shows  the  first 
several  columns  of  the  worksheet  “MATLAB  Results  for  Graphs”  which 
holds  the  meshed  raw  data  (Traw),  as  well  as  the  cleaned-up  raw  data 
(Tcleaned),  cleaned-up  raw  data  with  estimated  missing  values  (Tfilled), 
and  filtered  data  using  cutoff  periods  of  1,  2,  and  4  hours  (Tfiltl,  Tfilt2,  and 
Tfilt4,  respectively). 
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Figure  7b.  Example  output  from  CIRCAD.  This  figure  shows  a  results 
column  in  the  worksheet  named  “MATLAB  Results  for  Graphs”.  This 
column  is  the  result  of  fitting  cosine  curves  to  36-hour  data  windows, 
allowing  the  cosinor  function  to  fit  the  period,  and  working  with  unfiltered 
data.  Results  for  times  13.92  to  23.98  are  set  to  zero  (the  missing  value 
code)  because  curve  fitting  starts  at  midnight  of  the  first  day  of  data 
collection. 
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Figure  8.  Example  output  from  CIRCAD.  This  is  a  portion  of  the  “MATLAB 
Summary  Results”  worksheet  showing  the  fit  parameters  that  result  from 
fitting  cosine  curves  to  36-hour  data  windows,  allowing  the  cosinor 
function  to  fit  the  period,  and  working  with  unfiltered  data. 
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APPENDIX  A 

FILE  LIST 


Circad.m 

CircadianDataPath.m 

CosinorTpred.m 

DecodeCircadFileName.m 

DecodeJulian.m 

DoCosinor.m 

FFT_Filtering.m 

GetFileHandle.m 

InitializeDDE.m 

Lomb.M 

LombPeriod.m 

MeshData.m 

PokeTimeBasedColumn.m 
ReadCircadData.m 
RemoveOutliers.m 
Sumsqerror.m 
SumsqerrorFP.m 
T  estDDEConnection.m 
T  ransformSIeepWake.m 
WriteSummaryResultsByMethod.m 
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APPENDIX  B 

CIRCAD.M 


% - - - - - 

%Circad.m  ver  1.12  3  DEC  1999 

% - 

function  circad 

%  CIRCAD  limited  version  1.12  by  Tammy  Doherty,  USARIEM 
%  Analyzes  circadian  pill  temperature  data. 

%get  rid  of  old  variables  and  graphs 
ClearWorkspace 

%Read  data  from  .mat  or  excel  file 
X=ReadCircadData ; 

%Mesh  BCTM,  Minimitter  and  sleep/wake  data  to  common  date/time  reference 
X=MeshData (X) ; 

%Transform  Sleep/Wake  data  to  numeric  values  (38.2  for  sleep,  0  for  awake) 
X=Trans forms 1 eepWake (X) ; 

%Remove  outliers  from  pill  temperature  data 
X.Tclean=RemoveOutliers (X.Traw) ; 

%  Filter  Data  to  exclude  periods  less  than  1,  2,  or  4  hours. 
X=FFT_Filtering (X) ; 

%Specify  analyses  to  run: 

%In  this  version,  only  one  Analysis  type  is  available: 

AnalysisType (1) ={ 'cos inor  fit'}; 

Analysis Function (1) ={ 'DoCosinor ' } ; 

%Perform  analyses  using  raw  or  filtered  data, 

%using  data  intervals  of  various  widths  and  starting  locations, 

%and  using  various  methods  of  obtaining  circadian  period. 

RunAnalyses (X, AnalysisType, AnalysisFunction) ; 

return 

% - 

% - FUNCTION  DEFINITIONS - 

% - 


function  ClearWorkspace 
close  all; 
clear  all; 
return 

% - 

function  RunAnalyses  (X,  AnalysisType,  AnalysisFunction) 
for  i=l :prod(size (AnalysisType) ) 

disp (strcat (' starting  AnalysisType, '  ...')); 
disp ( '  ' ) ; 

[Var, VarDesc] =GetVar2Use (X) ; 
while  length (Var) >0 

[PeriodStr, PeriodDesc) =GetPeriodType2Use (X)  ; 
while  length (PeriodStr) >0 

Results^ [ ] ; 

[  Interval ,  IntervalDesc ,  Day,  RelDay ,  LastDayYN,  AppendYN]  =GetInterval2Use  (X) 
Per iod=AssignNumericCodeTo Period ( PeriodStr , X, Var , Interval ) 

while  length (Interval .XI) >0 

Title=strcat (X. FileHandle, ' :  ', char (AnalysisType ( i ) ), '  to',... 

VarDesc , '  over ' , IntervalDesc , PeriodDesc ) ; 
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%get  the  results  for  the  specified  interval 
IntervalResults=feval (char (AnalysisFunct ion (i) ) ,  .  .  . 

X.Time (Interval .Analysisl : Interval .Analysis2 ) .... 

Var (Interval .Analysisl : Interval. Analysis2) , Period) ; 

%if  results  need  to  be  appended  together,  do  it. 

%Otherwise,  the  final  results  array  equals  the  current  array 
if  AppendYN==l 

Results^AppendResults  (Results,  IntervalResults  ,X,  Day,RelDay,  Interval)  ; 
else 

Results=IntervalResults ; 

end 

%new  fields  (added  outside  of  the  analytical  routines) 
if  (AppendYN==l&Day==l) |AppendYN==0 

Results . StartDayRelative2Event=RelDay; 

Results .Title=Title; 

Results .t=X. Time; 

end 

%put  spacers  at  the  beginning  of  the  appended  array  so  the  length 
%and  times  correspond  to  X.Time 
if  -isempty (Day) 

if  Day==l  . . . . -  - - -  — ■  -  • 

%fill  all  arrays  up  to  il  with  NaN 
Results=InitResults (Results, Interval .Appendl-1) ; 

end 

end 

%Output  results 

if  (AppendYN==0) | (LastDayYN==l) 

%save  summary  results  to  an  Excel  worksheet 
WriteSummaryResultsByMethod (Results ) ; 

%plot  results 

PlotResults (Results, X, Var, VarDesc) ; 

%print  graph  to  file 

GraphFi leName=s treat (X.FileHandle,  '_Period# , num2str  (Period)  ,  ,  IntervalDesc,  ,  VarDesc, 

' .jpg' ) ? 

print ( ' -djpeg90 ' , fullf ile (X. Pa thName, GraphFi leName) ) ; 

%save  time-based  results  to  an  Excel  worksheet 
PokeTimeBasedColumn (Results .Title, Results . Tpred) ? 

end 


[Interval,  IntervalDesc,  Day, RelDay,  LastDayYN,  AppendYN]  =GetInterval2Use  (X) ; 
end 

[PeriodStr, PeriodDesc] =GetPeriodType2Use (X) ; 

end 

[Var, VarDesc]  =GetVar2Use (X)  ; 
end 

end 

disp ( ' DONE ! ') ; 
return 

% - - - 

function  PlotResults (Results, X, Var, VarDesc) 
persistent  FigNum 

if  is empty (FigNum) 

FigNum- 3 ; 

end 

N=length (Results. Tpred) ; 
i=find (Results .Tpred==0) ; 

Results .Tpred (i) =NaN; 
figure (FigNum) ; 

subplot  (2, 1, 1)  , plot  (X.Time, X.T,  'b' )  ,ylabel  ( 'unfiltered  Tpill' )  ,  title  (Results  .Title)  ; 
hold  on 

%plot  (X.Time,  Var,  'k' )  ; 

plot (X.Time (1 :N) , Results .Tpred (1 :N) , #r . ' ) 
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hold  off 

subplot (2, 1,2) /plot (X.Time, Var, 'b' ) ,ylabel (VarDesc) ; 
hold  on 

plot (X.Time (1:N) , Results . Tpred ( 1 :N) , 'r. ' ) 
hold  off 

%FigNum=FigNum+l ; 
return 


% - 

function  [composite] =AppendRe suits (composite, OneDay, X, Day, RelDay, Interval) 
if  Day=-1 
composite= [ ] ; 
end 

composite . t ( Interval . XI : Interval . X2 , 1 ) = . . . 

X . Time ( Interval . XI : Interval . X2 ) ; 
composite . Tpred ( Interval . XI : Interval . X2 , 1 ) = . . . 

OneDay . Tpred ( Interval . Appendl : Interval . Append2 ) ; 
composite.R(Day, 1) =OneDay.R; 
compos  ite.N{  Day, l)=OneDay.N; 
compos ite.MaxData (Day, 1) =OneDay . MaxData ; 
compos ite.MinData (Day, 1) =OneDay .MinData; 
composite.MaxPred(Day,  1)  =OneDay.MaxPred; 
compos ite.MinPred( Day,  1) =OneDay .MinPred; 

composite .TimeOf MaxData (Day, 1) =OneDay. TimeOf MaxData;  . 

compos  ite.TimeOfMinData  (Day,  1)  =OneDay.  TimeOf  MinData; 
compos  ite.TimeOfMaxPred  (Day,  1)  =OneDay.TimeOfMaxPred; 
composite.  TimeOfMinPred  (Day,  1)  =OneDay. TimeOf MinPred; 
compos ite.mesor (Day, 1) =OneDay.mesor ; 
composite .  amplitude  (Day,  1)  =OneDay  .amplitude; 
composite. period (Day, 1) =OneDay.period; 
composite. phase (Day, 1) =OneDay .phase; 
composite .gphase (Day, l)=OneDay.gphase; 
compos  ite.clock_j>hase  (Day,  1)  =OneDay.clock_phase; 
composite. phase_clockstring (Day, : ) =OneDay.phase_clockstring ( : ) ; 
composite. DayRelative2 Event (Day, 1) =RelDay; 
return 

% - 

function  [composite] =InitResults (composite, il) 
composite.Tpredd:  il,  1)  =NaN; 
composite. t (1; il, 1) =NaN; 
return 

% - 

function  [Var, VarDesc] =GetVar2Use (X) 
persistent  VarNum 
NumVars=4 ; 
if  is empty (VarNum) 

VarNum=l; 

else 

VarNum=VarNum+l ; 
end 

if  VarNum<=NumVars 
switch  VarNum 
case  1 

Var=X.Tclean; 

VarDesc=/Tpill_cleaned' ; 
case  2 

Var=X. Tf iltl ; 

VarDesc='  Tpill_f ilt (p=lhr) 7 ; 
case  3 

Var=X.Tfilt2; 

VarDesc='  Tpill_f ilt (p=2hr) ' ; 
case  4 

Var=X.Tfilt4; 

VarDesc=/  Tpill_f ilt (p=4hr) • ; 
case  5 

Var=X.Twrist; 

VarDesc='  Twrist__raw'  ; 
end 
else 

VarNum = [ ] ; 

Var=  [  ]  ; 
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VarDesc= ( ] ; 
end 
return 

% - 

function  [ Interval ,  IntervalDesc,  Day,  RelDay,  Las tDayYN,  AppendYN]  =GetInterval2Use (X) 

persistent  IntervalNum 

persistent  Numlntervals 

persistent  Circadianlnterval 

persistent  FirstMidnightlndex 

persistent  FirstNoonlndex 

persistent  NumM2MDays 

persistent  NumN2NDays 

persistent  N 

if  isempty (IntervalNum) 

IntervalNum=l; 

N= 1 eng th(X. Time) ? 

CircadianInterval=round(24/X.dTime) ; 

FirstMidnightIndex=round( (24-X.Time(l) ) /X.dTime+1) ; 
while  FirstMidnightIndex<l 

FirstMidnightIndex=FirstMidnightIndex+CircadianInterval; 

end 

LastMidnightIndex=FirstMidnightIndex+CircadianInterval; 
while  Las tMidnightIndex<=N 

LastMidnightIndex=LastMidnightIndex+CircadianInterval; 

end 

LastMidnightIndex=LastMidnightIndex-CircadianInterval; 

FirstNoonIndex=round(  ( 12 -X. Time (1)  )/X.dTime)  ; 
while  FirstNoonlndexcl 

FirstNoonIndex=FirstNoonIndex+CircadianInterval; 

end 

LastNoonIndex=FirstNoonIndex+CircadianInterval ; 
while  LastNoonIndex<=N 

LastNoonIndex=LastNoonIndex+CircadianInterval; 

end 

LastNoonIndex=LastNoonIndex-CircadianInterval ; 

NumM2MDays=f loor ( (LastMidnightlndex-FirstMidnightlndex+l ) /Circadianlnterval) ; 
NumN2NDays= floor ( (LastNoonlndex-FirstNoonlndex+l) /Circadianlnterval) ; 
Numlntervals  =NumM2MDays  *2  +NumN2NDays+l  ; 
else 

IntervalNum=IntervalNum+l  ; 

end 

LastDayYN=0 ; 

if  IntervalNum<=NumIntervals 

if  IntervalNum==l  %Do  Analysis  and  Lomb  on  whole  curve 
Day= [ 3 ; 

RelDay= [ ] ; 

AppendYN=0; 

Interval. Xl=l; 

Interval .X2=N; 

Interval .Analysis 1=1; 

Interval . Analysis2=N; 

Interval . Lombl=l  ; 

Interval . Lomb2=N; 

Interval.Appendl=l; 

Interval . Append2=N; 

IntervalDesc='  all'; 

elseif  (IntervalNum>l) & (IntervalNum<=NumM2MDays+l) 

%Do  Analysis  M2M,  Lomb  +/-12  from  Midnight-to-Midnight 
Day=IntervalNum-l ; 
if  Day = =NuinM2 MDay s 
Last Day YN=1  ; 

end 

AppendYN=l; 

RelDay=Day+X.StartDateNum- floor (X. Event Da teNum) ; 

Interval .XI =FirstMidnightIndex+ (Day- 1) *CircadianInterval; 

Interval . X2=FirstMidnightIndex+Day*CircadianInterval ; 

Interval .Analysisl=Interval.Xl; 

Interval . Analysis2 =Interval ,X2 ; 
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Interval .Appendl=l; 

Interval .Append2=Interval .X2- Interval .Xl+1; 

[Interval.Lombl,  Interval .Lomb2 ]  =GetOff setlnterval  — 

( Interval . XI , Interval . X2 , Circadianlnterval , N, 12) ; 

IntervalDesc='  24-hour  (midnight- to-midnight) ' ; 
el  seif  (IntervalNum>l+NumM2MDays)  &  (IntervalNum<=NumM2MDays+NumN2NDays+l) 
%Do  Analysis  N2N,  Lomb  +  /-12  from  Noon-to-Noon 
Day=IntervalNum-NumM2MDays-l ; 
if  Day==NumN2NDays 
LastDayYN=l; 

end 

AppendYN=l ; 

Re 1 Day=Day+X. S tar tDa teNum- floor  (X.  Event  Da  teNum)  ? 

Interval. Xl=FirstNoonIndex+  (Day-1)  *CircadianInterval; 

Interval .  X2=FirstNoonIndex+Day*CircadianInterval  ; 

Interval . Analysisl=Interval .XI; 

Interval  .Ana lysis2=Interval  ,X2  ; 

Interval  .Appendix- 

Interval  . Append2  =Interval . X2 -Interval . Xl+1 ; 

[  Interval .  Lombl ,  Interval .  Lomb2  ]  =GetOf  f  setlnterval  — 

( Interval . XI , Interval . X2 , Circadianlnterval , N,  12 )  ; 

IntervalDesc='  24-hour  (noon-to-noon) '  ; 
elseif  ( IntervalNum>l+NumM2MDays+NumN2NDays ) 

%Do  Analysis  Spin-to-Gam,  Lomb  +/-12  from  Midnight-to-Midnight 
Day=  IntervalNum-NumM2MDays  -NumN2NDays  - 1  ; 
if  Day==NumM2MDays 
LastDayYN=l; 

end 

AppendYN=l; 

RelDay=Day+X.  StartDateNum- floor  (X.  Event  Da  teNum)  ; 

Interval. Xl=FirstMidnightIndex+  (Day-1)  ^Circadianlnterval ; 

Interval  .X2=FirstMidnightIndex+Day*CircadianInterval; 

[ Interval. Analysisi,  Interval. Analysis2]=Get0ff setlnterval.  .  . 

( Interval . XI , Interval . X2 , Circadianlnterval ,  N,  6 ) ; 

Interval  .Appendl= Interval  .Xl-Interval  .Analysisl+1; 

Interval .  Append2= Interval  .Appendl+ Interval .  X2- Interval  .XI  ; 

[ Interval. Lombl, Interval .Lomb2] =GetOffse tlnterval . . . 

( Interval . XI , Interval . X2 , Circadianlnterval ,  N,  12 ) ; 

IntervalDesc='  36-hour  (6pm-to-6am) ' ; 

end 

else 

IntervalNum= [ ] ; 

Interval .Xl= [ ) ; 

Interval .X2= [ ) ; 

IntervalDesc=[] ; 

Day= [ ] ; 

RelDay= [ ] ? 

LastDayYN= [ ] ; 

AppendYN= [ ] ; 
end 
return 

% - 

function  [Period, PeriodDesc] =GetPeriodType2Use (X) 
persistent  PeriodType 
NumPeriodTypes=3  ; 
if  isempty (PeriodType) 

PeriodType=l; 

else 

Per iodType=PeriodType+l ; 

end 

if  PeriodType<=NumPeriodTypes 
switch  PeriodType 
case  1 

Period=/ free' ; 

PeriodDesc=' 
case  2 

Period= ' 24hours ' ; 

PeriodDesc='  period  set  to  24  hours'; 
case  3 

Period='Lomb' 

PeriodDesc='  period  set  to  Lomb  period'; 
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end 

else 

Period= [ ] ; 

PeriodType= [ ] ; 

PeriodDesc= [ ]  ; 

end 

return 

% - 

function  [Of fsetlndexl,0ffsetlndex2] =GetOf fsetlnterval. . . 

{ Interval_Xl ,  Interval_X2 ,  Circadianlnterval ,  N,  Of  f  setHours) 

fOf f set=Of fsetHours/24; 

Of fsetIndexl=Interval_Xl-fOf f set*CircadianInterval; 
if  OffsetIndexl<l 
Of fsetlndexl=l; 

end 

Of f setIndex2=Interval_X2+fOf f set*CircadianInterval; 
if  Of f setIndex2>N 
0ffsetIndex2=N; 

end 

return 

% - 

function  [Period] =AssignNumericCodeToPeriod ( Per iodStr#X, Var , Interval) 
if  strcmp (PeriodStr, ' free' ) 

Period=0; 

elseif  strcmp (PeriodStr, /24hours') 

Period=24; 

elseif  strcmp (PeriodStr, 'Lomb' ) 

Per iod=LombPer iod  (X . Time  ( Interval .  Lombl :  Interval .  Lomb2 )  -X . Time  ( Interval .  Lombl ) 
Var ( Interval . Lombl : Interval . Lomb2 ) ) ? 

end 

return 
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APPENDIX  C 

READCIRCADDATA.M 


% - 

%ReadCircadData.m  3  DEC  1999 

% - 

function  [X] =ReadCircadData 
%ReadCircadData.m  by  Tammy  Doherty,  USARIEM 
%reads  circadian  data  from  an  Excel  file 

%assumes  the  header  row  containing  variable  names  will  be  within  the  first  10  rows  and  30 
columns 

%get  file  name 

[X . PathName , X . FileName , X . FileHandle] =GetFileName; 

%decide  between  the  excel  file  and  mat  file  (if  it  exists) 

[ReadMatFileYN]  =DecideBetweenExcelAndMatFile  (X. PathName, X. FileHandle)  ; 

if  ReadMatFileYN==l 
X=ReadMatData (X) ; 
else 

X=ReadExcelData (X) ; 

%save  data  in  a  binary  '.mat'  file  for  next  time 
Save2MatFile (X) ; 

end 

return 


function  [ReadMatFileYN] =DecideBetweenExcelAndMatFile (PathName, FileHandle) 
MatFileName=s treat (FileHandle, ' .mat' ) ; 

ExcelFileName=strcat (FileHandle, ' .xls' ) ; 

DirList=dir (PathName) ; 

MatFileListNum=0;  ExcelFileListNum=0; 

for  i=l : length (DirList) 

if  stremp (DirList (i) .name,MatFileName) 

MatFileListNum=i ; 

el seif  stremp (DirList (i) .name, ExcelFileName) 

ExcelFileListNum=i ; 

end 

end 

if  MatFileListNum>0 
ReadMatFileYN=l ; 
else 

ReadMatFileYN=0 ; 
end 

return 

% - 

function  [X] =ReadMatData (X) 

TestDDEConnections (X. FileHandle) ; 

MatFileName=strcat (X. PathName, X. FileHandle) ; 

X=load(MatFileName) ? 
return 

% - 

function  [X] =ReadExcelData (X) 

%notify  user  that  data  read  is  starting 

disp( 'reading  data  (this  could  take  several  minutes)  ...'); 
disp  ( '  '  )  ; 

TestDDEConnections (X. FileHandle) ; 

%Read  BCTM,  MiniMitter,  and  Sleep /Wake  data 
X=ReadDeviceData (X, 'BCTM' , 'PillTemp' , 'Pill  Temp' ,1,1,1); 
X=ReadDeviceData (X, 'Mini' , [ 'Twrist  ' ; 'Tother  ' ; 'Activity' ) , _ 
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['Twrist  ';'Tother  '; 'Activity' ], [1  1  1] , 0, length (X.BCTM„Hour) ) ; 
if  length (X.Mini_Hour) ==0 
X . Mini_Hour=X . BCTM_Hour ; 

end 

X=ReadDeviceData (X, 'Sleep' , 'SleepWake' # 'Sleep/Wake' ,0,0, length (X.BCTM_Hour) ) ; 
if  length (X. Sleep_Hour) ==0 
X . Sleep_Hour=X . BCTM_Hour ; 

end 

%Read  General  Information  (i.e.,  LH  peak.  Cycle  Start,  Data  Start  and  Data  Stop) 
X=ReadGenInf o (X) ? 

return 

% - 

function  TestDDEConnections (FileHandle) 

%TestDDEConnection  is  located  in  TestDDEConnection.m  and  is  responsible 
%for  displaying  an  error  box  and  stopping  program  execution  if  the  "topic" 

%is  not  found  {i.e.,  assigned  channel  number  is  <=0) 

%Test  the  raw  data  worksheet  is  there  (same  name  as  the  file  name) 
TestDDEConnection (FileHandle) 

%test  whether  the  'Genlnfo'  worksheet  is  there 
TestDDEConnection ( 'Genlnfo' ) 

%test  whether  the  'MATLAB  Summary  Results'  worksheet  is  there 
TestDDEConnection ( 'MATLAB  Summary  Results') 

%test  whether  the  'MATLAB  Results  for  Graphs'  worksheet  is  there 
TestDDEConnection ( 'MATLAB  Results  for  Graphs') 
return 

% - - 

function  CloseDDEConnect ion (channel) ; 

dde term (channel) ; 
return 

% - 

function  [channel] =OpenDDEConnect ion (FileHandle) ; 

channels Ini tializeDDE (FileHandle) ; 
return 

% - - - 

function  [X]  =ReadDeviceData  (X,  Device,  FieldName,  ColHeading,  VarType,RequiredYN,  AltN)  ; 
channel=OpenDDEConnection (X . FileHandle) ; 
if  length (Device) >0 

Time Var=s treat (Device, '_Time' ) ; 

DayVar=s treat (Device, '_Day' ) ; 

DateVar=s treat (Device, '_Date' ) ; 

ElDaysVar=strcat (Device, '_ElDays' ) ; 

TimeColHead=s treat (Device, '-Time' ) ; 

DayColHead=s treat (Device, '-Day' ) ? 

DateColHead=strcat (Device, ' -Date' ) ; 

ElDaysColHead=s treat (Device, ' -ElDays' ) ; 
else 

TimeVar= ' Time ' ; 

DayVars'Day'  ; 

TimeColHead='Time' ; 

DayColHead=='Day'  ; 

end 

eval (streat { 'Time=ReadExcelColumn( channel, ' ' ' ,TimeColHead, ' ' ' , 50000, 0,0);')) 
if  length (Time) ==0&RequiredYN==l 

errordlg (streat ( 'Must  Enter  Either  Time  or  ',TimeColHead),'ERROR'); 
error (streat ( 'Must  Enter  Either  Time  or  ' , TimeColHead) ) ; 
elseif  length (Time ) ==0&RequiredYNs=0 
N=AltN; 

elseif  length (Time ) >0 
[N, dummy] =size (Time) ; 

end 

eval (streat ( ' Day =ReadExcelColumn (channel , ' ' ' , DayColHead, ' ' ' ,N, 1, 0) ; ' ) ) ; 
if  ~isempty{DateColHead) 

eval (streat ( ' Date=ReadExcelColumn (channel , ' ' ' , DateColHead, ' ' ' ,N, 0, 0) ; ' ) ) ? 

end 
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if  -isempty (ElDaysColHead) 

eval (s treat ( ' ElDays =ReadExcelColumn (channel, ' ' ' , ElDaysColHead, ' ' ' ,N,  1,  0) ; ' ) )  ? 

end 

if  length (ElDays) ==0&length{Day) ==0&length(Date) ==0&RequiredYN==l 

errordlg (s treat ( 'Must  Enter  Day,  ' , DayVar, # ,  ' , DateVar , ' ,  or  ' , ElDaysVar ) , ' ERROR' )  ; 
error ( s treat ( 'Must  Enter  Day,  DayVar,',  ', DateVar,',  or  ', ElDaysVar) ) ; 

end 

[NumFields ,  dummy]  =size  (FieldName)  ; 
for  i=l : NumFields 

eval  { s treat  ( 'X.  ' ,  FieldName  (i,  : ) ,  '  =ReadExcelColumn( channel,  ' '  ' ,  ColHeading  (i,  :),... 

"  '  ,N,  ' , num2str (VarType { i ) ) , ' ,0) ; ' ) ) ; 

end 

CloseDDEConnection (channel) ; 

disp  (s  treat  ( 'Done  ^Device,'  read...')); 

%convert  Time/Date  information  to  'Device_Datenum'  and  'Device_Hour '  fields 

eval  (s  treat  ('X.  '#Device,  '  _DateNum=Convert  2  Da  teNum(  Day,  Date,  ElDays,  Time,  X.FileName)  ;')); 
eval  (s  treat  ( 'X. Device,  '_Hour=Convert2Hour  (X. Device,  '_DateNum,  X.  FileName)  ;')); 
disp (s treat ( 'Done  ', Device,'  Date/Time  conversions ...')) ; 

return 

% - - 

function  [X] =ReadGenInfo (X) 

[channel] =InitializeDDE ( 'Genlnfo' ) ; 

X.StartDateNum=GetStartDateNum (X.FileName) ; 

X.  Event  Da  teNum=datenum(ddereq(  channel,  'r5c2' ,  [1  1] )  )  ? 

Da taS tartDateString= (ddereq( channel, 'r7c2' , [1  1] ) )  ; 

if  double  (DataStartDateString)  ==10  %condition  upon  reading  a  blank  cell 
X. DataStartDateNum=X . StartDateNum; 
else 

X . DataStartDateNum=datenum (DataStartDateString ) ? 

end 

DataStopDateString= (ddereq( channel, 'r8c2 ' , [1  1] ) ) ; 

if  double (DataStopDateString) ==10  %condition  upon  reading  a  blank  cell 
%use  the  date  information  from  the  last  entry  from  the  BCTM  data 
X . DataStopDateNum=f loor (X . BCTM_DateNum  ( length (X . BCTM_DateNum) ) ) ; 

else 

X . Da taStopDateNum=datenum (DataStopDateString)  ; 

end 

ddeterm( channel) ; 
return 


function  [PathName, FileName, FileHandle]  =GetFileName 

%get  the  current  directory  name 
working__directory=cd; 

%set  current  directory  to  Circadian  Data  Path  (from  CircadianDataPath.m) 
cd(CircadianDataPath) ; 

%get  data  file  name 

[FileName,  PathName]  =uigetf  ile  ( '  *  .xls ' ,  'select  data  file'); 

FileHandle=GetFileHandle (FileName)  ; 

%change  back  to  original  working  directory 
cd(working_directory)  ; 

return 

% - 

function  (Row,  Col]  =Get Star tRowCol  (channel,Vname) 

%Assumes  the  header  row  containing  variable  names  will  be  within  the  first  10  rows 
%and  30  columns 
Row=0;  Col=0; 
maxrows=10;  maxcols=30; 
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%find  the  row  and  column  numbers  of  the  recognized  variable  names 
row=l ,- FoundMatch=0 ; 
while  row<=maxrows&~FoundMatch 
col=l; 

while  col<=maxcols&~FoundMatch 

x=ddereq (channel, s treat ( 'r' , num2str (row) , 'c' , num2str (col) ) , [1  1]  )  ; 
x=x (1: length (x) -1) ; 
if  stremp (upper (x) , upper (Vname) ) ==1 
Row=row;  Col=col; 

FoundMatch=l; 

end 

col=col+l; 

end 

row=row+l; 

end 

return 


function  [X]  =ReadExcel  Column  (channel,  Vname ,  N,  Nume  r  i  cYN,  checkN) 

%string  format  (sformat)  for  reading  from  Excel 
% Indx=Vname Index ( Vname ) ; 

[startrow, col] =GetStartRowCol (channel, Vname) ; 
if  (startrow~=0) & (col~=0) 
lastrow=startrow+N; 

Xcel larray=ddereq( channel, s treat ( 'r' , num2str (startrow+1) , 'c' ,num2str (col) , . 

' :r' ,num2str (las trow) , 'c' ,num2str (col) ) , [1  1]); 

Xstringarray=char (Xcellarray) ; 
i=findstr (Xstringarray, char (10) ) ; 

N=length ( i ) ; 
il=i (1 :N-1) +1; 
il= [1  il]; 
i2=i-l ; 
if  NumericYN 

X ( 1 :N) =NaN;  %initialize  to  NaN  in  case  of  missing  values 
X=X' ; 

end 

for  z=l:N 

if  il (z) <=i2 (z) 

Xstringmatrix ( z , 1 : (i2 (z) -il (z) +1) ) =Xstringarray (il (z) : i2 (z) ) ; 
if  NumericYN& length (Xstringmatrix ( z , :))>0 
if  length ( str2num (Xstringmatrix ( z , : ) ) ) >0 
X (z) =str2num (Xstringmatrix (z, : ) ) ; 

end 

end 

end 

end 

if  -NumericYN 

if  isempty (Xstringmatrix) 

X=  [  ]  ; 

N=0  ; 
else 

X=Xstringmatrix; 

end 

end 

else  %case  if  the  variable  is  not  present  in  the  file  (row=0,  col=0) 

X=  ( ]  ; 
end 

%pad  end  of  array  to  get  N  elements 
if  checkN==l; 

[rows, cols] =size (X) ; 
if  rows>0 
for  row=rows+l:N 
if  NumericYN 

X (row, : ) =NaN ; 
else 

X (row,  : ) ={ ' ' } ; 

end 

end 

end 

end 

return 
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function  [DateNum] =Convert 2 DateNum (Day, Date, ElapsedDays , Time, fname) 

[ Sub j  ect, MPhase , FileNum, MCycle , Julian, St art Year ] =DecodeCircadFileName ( fname) ; 
[StartMonth, StartDay] =Decode Julian (Julian, StartYear) ? 
if  length (Date) >1  %date  entered  as  string 
DateNum=datenum(Date)  +datenum(Time)  ; 
else 

if  length (Day) >1 
N=  length (Day) ; 

Year ( 1 : N, 1 ) =S tartYear ; 

Month (1:N, 1) =StartMonth; 

iMonthChange=f  ind ( (datenum(Year  (2  :N)  ,Month(2:N)  , Day(2:N) ) -datenum(Year  (1  :N- 

1),... 

Month (1 :N-1) , Day (1 :N-1) ) ) <0) ; 
while  length ( iMonthChange ) >0 

NewYear=Year ( iMonthChange ( 1 ) ) ; 

NewMonth=Month ( iMonthChange ( 1 ) ) +1 ; 
if  NewMonth>12 
NewMonth=l; 

NewYear=NewYear+l ; 
end 

Month ( iMonthChange ( 1 )  + 1 :  N ) =NewMonth ; 

Year ( iMonthChange (1) +1 :N) =NewYear; 

iMonthChange = f ind ( (datenum(Year (2 :N) ,Month(2:N) ,Day(2:N) ) -datenum(Year (1 ;N- 

1),... 

Month ( 1 : N- 1 ) , Day ( 1 : N- 1 ) ) )<0) ; 

end 

DateNum=datenum (Year, Month, Day)  +datenum(Time) ; 
elseif  length (ElapsedDays) >1 

DateNum=datenum (StartYear , Star tMonth , StartDay) +ElapsedDays+datenum(Time) ; 
else 

DateNum= [ ) ; 

end 

end 

return 


function  [Hour] =Convert2Hour (DateNum, fname) 
if  isempty (DateNum) 

Hour= [ ] ; 
else 

StartDateNum=GetStartDateNum ( fname) ; 

Hour= (DateNum-StartDateNum) *24; 

end 

return 

% - 

function  [StartDateNum] =GetStartDateNum( fname) 

[Subject, MPhase, FileNum, MCycle Julian, StartYear] =DecodeCircadFileName ( fname) ; 
[StartMonth, StartDay] =DecodeJulian (Julian, StartYear) ; 

StartDateNum=datenum (StartYear, StartMonth, StartDay) ; 

return 

% - 

function  Save2MatFile (X) 

MatFileName=strcat (X . PathName, X. FileHandle) ; 

FieldName=f ieldnames (X) ; 

FileDate=datestr (now, 0); 
save (MatFileName, 'FileDate' ) ; 
for  i=l: length (FieldName) 

eval (s treat (char (FieldName (i) ) , ' =X. # , char (FieldName (i) ),';')) ; 
save (MatFileName, char (FieldName (i) ) , '-append') ; 

end 

return 
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APPENDIX  D 

MESHDATA.M 


%MeshData.m  3  DEC  1999 

% - 

function  [Y] =MeshData (X) ; 

disp { 'meshing  variables  into  a  single  time  base 

%copy  non-meshed  fields  to  Y 
Y=CopyNonMeshedFields2Y(X) ; 

%Get  data  start  and  stop  times.  Start  times  are  either  the  first  time  encountered 
%in  any  of  the  device  date/time  fields,  or  the  data  start  time  as  entered  into  the 
%GenInfo  worksheet.  Stop  times  are  either  the  latest  time  encountered 
%in  any  of  the  device  date/time  fields,  or  the  data  stop  time  as  entered  into  the 
%GenInfo  worksheet. 

[Y. MinHour , Y . MaxHour ] =GetDataStartS topTimes (X)  ; 

%Compute  new  Hour  array 
Y.dTime=l/60;  % interval  is  1-minute 
Y.Time= [ Y. MinHour :Y.dTime :Y. MaxHour] ' ; 

Y =DoTheMapp ing ( X , Y ) ; 

%send  results  to  'MATLAB  Results  for  Graphs'  worksheet 
SendMeshedResultsToExcel (Y) ; 

return 

% - 

function  [Y] =CopyNonMeshedFields2Y  (X) 

MeshedFields='_DateNum,  _Hour ,  PillTemp ,  Twrist ,  Tother ,  Activity,  SleepWake '  ; 

FieldName= fieldnames  (X)  ; 
for  i=l : length (FieldName) 

if  isempty ( findstr {MeshedFields , char (FieldName ( i) )))... 

&isempty (finds tr( char (FieldName (i) ) , '_Hour' ))&... 
isempty (findstr (char (FieldName (i) ) , '_DateNum' ) ) 
eval (strcat ( ' Y. ' , char (FieldName ( i) ) , '=X. ', char (FieldName (i) 
end 

end 

return 

% - 

function  [MinHour, MaxHour} =GetDataStartStopTimes (X) 

StartHour= (X.DataS tartDateNum- floor (X.DataStartDateNum) ) *24; 

StopHour= (X. Da taStopDateNum- floor (X . DataStartDateNum) ) *24; 

%get  minimum  of  data  time  points 
MinHour =X . BCTM_Hour ( 1 ) ; 
if  -isempty (X.Mini__Hour) 

MinHour =min (MinHour , X .Mini_JH our  ( 1)  )  ; 
end 

if  -isempty (X.Sleep_Hour) 

MinHour =min (MinHour , X . Sleep_Hour ( 1 ) ) ; 
end 

%get  maximum  of  StartHour  and  MinHour 
MinHour=max (MinHour, StartHour) 

%get  maximum  of  data  time  points 
MaxHour =X . BCTM_Hour ( length (X . BCTM_Hour ) ) ; 
if  -isempty (X.Mini_H our) 

MaxHour =max( MaxHour, X.Mini_Hour( length (X.Mini_Hour) ) ) ; 

end 

if  -isempty (X.Sleep_Hour) 

MaxHour=max (MaxHour, X.Sleep_Hour (length (X. Sleep_Hour) ) ) ; 

end 

MaxHour =min ( S  topHour , MaxHour ) ; 

%round  MaxHour  and  MinHour  to  nearest  whole  minute 
MaxHour =round { MaxHour  *60)/60; 

MinHour =round (MinHour* 60) /60; 
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return 


function  [iStart, iStop] =GetStartStopElement (Device_Hour , StartHour, StopHour) 
i=find(Device_Hour<StartHour) ; 
if  isempty(i) 
iStart=l; 
else 

iStart=max { i ) +1 ; 
if  is  tar  t>  length  (Device__Hour) 
iS tar t= length (Device_Hour) ; 

end 

end 

i=find(Device_Hour>StopHour) ; 
if  isempty(i) 

iStop=length(Device_Hour) ; 
else 

iStop=min(i) -1; 
if  iStop<iStart 
iStop=iStart; 

end 

end 

return 

% - 

function  [X] =Round2NearestMinute (X) ; 

X . BCTM_RoundHour= (X . BCTM_Hour . *  60 ) / 6  0 
if  isempty (X.Mini_Hour) 

X . Mini_RoundHour = [ ] ; 
else 

X . Mini_RoundHour = (X.Mini__Hour .  *60) /60 

end 

if  isempty  (X.  Sleep__Hour) 

X . S 1 e ep_RoundHour = [ ] ; 
else 

X . Sleep_RoundHour= (X. Sleep_Hour . *60) /60 

end 

return 


% - 

function  [Y] =DoTheMapping (X, Y) 

%Compute  the  indices  for  mapping  the  old  arrays  into  the  new  arrays 
BCTM. i=round ( (round (X . BCTM_Hour . /Y. dTime) .*Y.dTime  -  Y.MinHour) ./Y.dTime)+l; 
Mini . i=round( (round(X.Mini__Hour . /Y.dTime) .*Y. dTime  -  Y.MinHour) ./Y.dTime)+l; 
Sleep. i=round( (round (X.Sleep_Hour./Y. dTime) .*Y. dTime  -  Y.MinHour) . /Y.dTime) +1 ; 

%Do  the  mapping 
N=length (Y.Time) ; 

%BCTM 

Y . Traw=Map (X . PillTemp , BCTM . i , N, 1 ) ; 

%MiniMitter 

Y.Twrist=Map (X.Twrist ,Mini . i,N, 1) ; 

Y . Tother=Map (X . Tother , Mini .i#N,l); 

Y. Activity=Map {X. Activity /Mini . i/N/ 1) ; 

%Sleep/Wake 

Y. SleepWake=Map (X. SleepWake, Sleep. i.N,  0) ; 
return 

% - 

function  [DataMapped]  =Map{Data,  iData.N^umericYN) 
if  length (Data) >0 

iGoodData=f ind( (iData>0) & (iData<=N) ) ; 

Data=Data (iGoodData) ; 
iData=iData (iGoodData) ; 
if  NumericYN==l 
DataMapped ( 1 :N) =NaN; 
end 

DataMapped ( iData ) =Data ;  DataMapped=DataMapped ' ; 
if  length (DataMapped) <N 
if  NumericYN==l 

DataMapped (N) =NaN; 
else 

DataMapped (N) =Data ( length (Data) ) ; 

end 
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end 

else 

DataMapped= [  ]  ; 

end 

return 

% - 

function  SendMeshedResultsToExcel (X) 

%poke  time  and  descriptive  variables  (pill#,  event#,...) 
PokeTimeBasedColumn ( ' Time ' , X . Time) ; 
PokeTimeBasedColumn ( 7 Traw7 , X . Traw) ; 
PokeTimeBasedColumn ( ' Activity' ,X. Activity) ; 
PokeTimeBasedColumn ( 7 Twr is t 7 ,X.Twrist) ; 
PokeTimeBasedColumn ( 7 Tother 7 ,X.Tother) ; 

return 
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APPENDIX  E 

DATA  TRANSFORMATION  FUNCTIONS 


% - 

%RemoveOutliers.m  10  AUG  1999 

% - 

function  [T] =RemoveOutliers (T) 

%set  outliers  to  NaN 

%note  that  index  will  also  identify  string  values 
%which  are  assigned  numeric  values  of  zero 
index=find( (T>40) | (T<35) ) ; 

% As sign  any  elements  of  T  that  are  in  index  the  value  NaN 
T  (index)  =NaN; 

%write  results  to  Excel  file 
PokeTimeBasedColumn( 'Tcleaned' ,T) ; 

return 


% - 

%TransformSleepWake.m  11  AUG  1999 

% - 

f unc t i on  [ X ] =Tr ans forms 1 eepWake ( X ) 

% transforms  SleepWake  from  a  string  array  to  a  numeric  array 

%where  "sleep"  is  given  a  value  of  38.2,  "wake"  is  given  a  value  of  zero 

%First,  look  to  see  whether  Sleep/Wake  was  already  converted 
NotConverted=isempty (f ind(X. SleepWake==38 . 2 ) ) ; 
if  NotConverted 

x=double (upper (X. SleepWake) ) ;  %converts  string  characters  to  ascii  code 
index=find(x==0) ; 
if  length ( index) >0 
if  index(l)==l 

x(l)=87;  %if  first  entry  is  blank,  set  to  "awake" 

end 

%for  any  blank  values,  drag  down  the  previous  value 
for  j =2 : length (index) 

x ( index ( j ) ) =x ( index ( j ) -1 ) ; 

end 

end 

%convert  all  code  83 's  to  38.2  and  87 ' s  to  zero 
IsAsleep= (x==83 ) ; 

x=IsAsleep. *38 .2;  %this  step  automatically  sets  non-83's  to  zero 
PokeTimeBasedColumn ( ' Sleep ' ,x) ; 

end 

return 
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APPENDIX  F 

FFT_FILTERING.M 


%FFT_Fi 1 tering . m  11  AUG  1999 

% - - 

% - 

function  [X] =FFT_Filtering (X) 

%f ill  in  blank  temperature  values  so  that  the  FFT  will  have 
%a  temperature  array  with  a  constant  sampling  interval 

%Reminder :  Tclean  is  the  meshed  pill  temperature  vector  with  text  and  outliers  removed 
X.T=FillTemperatureBlanks (X. Tclean) ; 

N= length (X.T) ; 

%N  needs  to  be  an  even  number 
if  (N/2-f loor (N/2) ) >0 
N=f loor (N/2) *2; 
end 

FFTofT=f f t (X.T (1  :N) ) ; 

%for  the  periodogram: 
freq={ (1/X.dTime) ,*(0:N/2) . /N) ' ; 
period (2 : floor (N/2 ) )  =  (1 . /freq(2 : floor (N/2 ) ) ) ' ; 
magnitude=abs (FFTofT) ; 

for  f_cutof f = [ . 25  .5  1.0] 

%set  to  zero  everything  with  freq>f_cutof f 
i_cutof f =ceil (f_cutof f / ( (1/X.dTime) ) *N) ; 

FFTf ilt=FFTofT; 

FFTf ilt (i_cutof f+1 :N/2 ) =0; 

FFTfilt (N/2+2 :N-i__cutof  f ) =0; 

%perform  inverse  fft  to  get  filtered  curve 
VarName=s treat ( 'X.Tf ilt' , num2str (round (l/f_cutof f ) ) ) ; 

Var=abs (if ft (FFTfilt) ) ; 
if  length (Var)< length (X.T) 

Var (N+l : length (X.T) ) =X.T(N+1: length (X.T) ) ; 
end 

%assign  the  filtered  curve  to  X.Tfiltl,  X.Tfilt2,  or  X.Tfilt4 
eval (streat (VarName, '=Var; ' ) ) ; 

end 


%plot  results 

FFT_plot ( freq, period, magnitude, X) ; 

%send  filtered  data  to  'MATLAB  Results  for  Graphs' 
PokeTimeBasedColumn( 'Tfiltl' /X.Tfiltl) ; 
PokeTimeBasedColumn ( 'Tf ilt2 ' , X.Tf ilt2 ) ; 
PokeTimeBasedColumn ( ,Tfilt4/ ,X.Tfilt4)  ; 

return 


function  FFT_plot ( freq, period, magnitude , X) 
figure (1) ; 

Np=length (period) ; 

%plot  raw  data 

subplot  (5, 1, 1)  ,  plot  (X. Time,  X.T)  ,ylabel  ( 'raw' )  ,xlabel  ( '  time' ) 
title (X.FileHandle) 

%plot  periodogram 

subplot (5, 1, 2) ,plot (period (2 :Np) , magnitude (2 :Np) ) ,ylabel ( 'magnitude' ) ,xlabel ( 'period' ) 
%plot  filtered  data 

subplot (5, 1,3) ,plot (X. Time, X.Tf iltl) ,ylabel (' filtered' ) , xlabel ( ' time' ) 
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Title=strcat (X. FileHandle, 7 :  period  cutoff  =  lhr'); 

subplot (5, 1,4) , plot (X.Time, X. Tf ilt2 ) ,ylabel { ' filtered' ) , xlabel ( ' time' ) 
Title=strcat (X. FileHandle, ' :  period  cutoff  =  2hr'); 

subplot (5, 1,5) ,plot (X.Time,X.Tfilt4) ,ylabel ( ' filtered' ) , xlabel ( ' time' ) 
Title=strcat (X. FileHandle, ' :  period  cutoff  =  4hr'); 
return 

% - 

function  [T] =FillTemperatureBlanks (T) ; 

NotANumber=isnan(T) ; 

IsZero= (T==0) ; 

IsBlank=NotANumber+IsZero; 
index=f ind(IsBlank>0) ; 
for  j=l: length (index) 
if  index(j)~=l 

T ( index ( j ) ) =T ( index ( j ) -1 ) ; 

end 

end 

% to  estimate  values  for  the  beginning  of  the  array,  go  backwards  from  1st  good  value 
NotANumber=isnan(T) ; 

IsZero= (T==0 ) ? 

IsBlank=NotANumber+IsZero; 
index=f ind(IsBlank>0) ; 
for  j =length( index) : -1 : 1 

T { index ( j ) ) =T ( index { j ) +1 ) ; 

end 

PokeTimeBasedColumn( 'Tfilled' ,T) 
return 
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APPENDIX  G 

FUNCTIONS  TO  GET  THE  LOMB  PERIOD 


% - 

%LombPeriod.m  3  DEC  1999 

% - 

%Name  changed  to  LombPeriod  from  LombAnalysis  on  1  SEP  1999 

function  [Period] =LombPeriod(t/ T) 
sig=0. 001; 

disp ( ' starting  Lomb  analysis  .  . . ' )  ; 
disp  ( '  ' )  ; 

%set  up  axes 

[xticks,xticklabels] =XAxesSetup (t) ; 

%Run  Lomb  on  Whole  Curve 

%The  lomb  function  used  here  was  written  by  John  Shimeld,  June  1997 
%obtained  at  www. swarthmore.edu/epigone/comp .  so  ft_„sys .matlab/shermgermsku/ 

%Pine.HPP. 3. 9 6. 971118082015. 305 6A-100000@agc2 
[ f Lomb, PerLomb, PowLomb, Psig] = lomb (t,T, [. 5  .05  .01  .001]); 

%plot  results 

goplot (PerLomb, PowLomb, 'period  (hours) ' , 'Lomb  Power' ,  'b' ,  2 , 1, 1, 1, xticks, xticklabels) ; 
goplot ( [ 0  max ( PerLomb) ] , [ Psig ( 1 )  Psig (1) ], 'period  (hours)',... 

'Lomb  Power' , 'y-' ,2, 1, 1, 1, xticks, xticklabels) ; 
goplot ([0  max (PerLomb) ], [Psig (2)  Psig (2) ], 'period  (hours)',... 

'Lomb  Power' , 'c- ' , 2, 1, 1, l,xticks ,xticklabels) ; 
goplot([0  max (PerLomb) ], [Psig (3)  Psig (3) ], 'period  (hours)',... 

'Lomb  Power' , 'g-' ,2, 1,1, l,xticks,xticklabels) ; 
goplot ( [0  max (PerLomb) ], [Psig (4)  Psig (4) ], 'period  (hours)',... 

'Lomb  Power' , 'r-' , 2, 1, 1, l,xticks, xticklabels) ; 
legend ( 'Lomb  Power' , 'P<0. 5 ' , 'P<0 . 05' , 'P<0 . 01' , 'P<0 .001' ) ; 

[MaxVal , Maxlndx] =max ( PowLomb) ; 

Peri od=PerLomb (Maxlndx) ; 
disp ('DONE  LOMB! ') ; 
return 


function  goplot  (xdata,ydata, xlabel, ylabel ,plotchar,  f  ignum,  numplothoriz,  .  .  . 
numplotvert,  plotnum, xticks, xticklabels) ; 
figure ( f ignum) ; 
grid  on; 
hold  on; 

subplot (numplotvert, numplothoriz, plotnum) ,plot (xdata,ydata,plotchar) ,  XLABEL (xlabel ) , 
YLABEL (ylabel ) ; 
hold  off; 

set (gca, 'XTick' , xticks) ; 

set(gca, 'XTickLabel' , xticklabels) ; 

return 

% - 

function  [xticks , xticklabels ] =XAxesSetup (t) 
maxtime= (floor (max (t) /24 ) + (mod(max (t) , 24) >0) ) *24; 
xticks= [0 :2 : maxtime] ; 
xticklabels=int2str ( 0 ) ; 
for  tick=6 : 6 :maxtime 

xticklabels=strcat (xticklabels, ' | | | ' ) ; 
xticklabels=strcat (xticklabels, int2str (tick) ) ; 

end 

return 
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%Lomb.m 

% - 

%  LOMB.M 
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%  Spectral  Analysis  of  Unevenly  Sampled  Data 

% 

%  Implements  the  Lomb  algorithm  (cf.,  Numerical  Recipes,  2nd  ed. , 

%  page  575)  to  estimate  the  spectral  power  of  an  unevenly 
%  sampled,  2D  dataset. 

% 

%  INPUT  PARAMETERS: 

%  t  -  sample  times  [nxl] 

%  z  measured  value  at  each  sample  time  [nxl] 

%  sig  -  significance  levels  for  which  a  power  estimate  is 

%  desired  (e.g.,  [0.5  0.05]  would  give  power  values 

%  for  the  50%  and  95%  confidence  limits) . 

% 

%  OUTPUT  PARAMETERS: 

%  f  frequencies  for  the  Lomb  power  estimates 

%  P  power  estimates 

%  Psig  -  powers  at  the  desired  confidence  limits 
% 

%  NOTE:  this  is  not  a  particularly  fast  implementation  of  the 
%  Lomb  periodogram.  See  Numerical  Recipes  for  a  faster 

%  algorithm. 

% 

%  John  Shimeld,  June  1997. 

% - - - 

function  [freq, period, power, Psig] =lomb(t, Z, sig) 

%remove  missing  values  from  Z  (modification  3  DEC  99  by  TJD) 
i=f ind(~isnan(Z) ) ; 
t-t (i) ; 

Z=Z ( i ) ; 

ofac  =  4;  %  oversampling  factor,  should  be  >=  4 

hifac  =  1/12;  %2;  %  determines  the  highest  frequency  for  which 

%  the  powers  are  estimated.  This  corresponds 
%  to  hifac* (Nyquist  frequency  if  dataset  were  evenly  spaced) 

N  =  length(t) ; 

Zmean  =  mean ( Z )  ; 

Zvar  =  cov ( Z )  ; 

t_span  =  max ( t ) -min ( t ) ; 

np  =  ofac*hifac*N/2 ; 

fc  =  N/2/t_span; 

fmin  =  l/t_span; 

fmax  =  hifac*fc; 

fstep  =  ( fmax- fmin ) /np; 

nf  =  ( fmax- fmin) /fstep+1; 

freq  =  [fmin: fstep: fmax] ; 
pstep= (max(t) -min(t) ) /np; 
period=min (t) :pstep:max ( t) ; 

for  z=l:np+l 

%This  modification  added  8/10/99  by  TJD  to  avoid  div-by-zero  error  in  omega  assignment 
statement 

if  period(z)==0 

period (z)=0. 0000000001; 
end 

omega  =  2 . *pi/period{z) ; 
tau  =  atan (sum (sin (2 . *omega. *t) )  .  /  ... 
sum(cos (2 . * omega. *t) ) ) . / 2 . /omega; 
a  =  cos (omega .* (t- tau) ) ; 
b  =  sin (omega. * (t-tau) ) ; 

power (z)  =  (sum( (Z-Zmean) . *a) . A2 . /sum(a. A2)  +  sum( (Z- 
Zmean) . *b) .^2 . /sum(b. A2) ) . / (2*Zvar) ; 
end 

Psig  =  log (hifac . *np. /ofac . /sig) ; 
return 
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APPENDIX  H 

DOCOSINOR.M 


%DoCosinor .m  4  OCT  1999 


function  [cosinor] =DoCosinor (t, T, Period) 

%need  to  remove  missing  values 
i=find(~isnan(T) ) ; 

T2=T{i) ; 
t2=t (i) ; 

if  Period==0 

p0= [40  .5  1  .25] ; 

p=fmins { ' sumsqerror ' , pO , [0] , [ ] ,T2, t2) ; 
elseif  Period==24 
p0= [40  .5  13; 

p=fmins { ' sumsqerrorFP ' , pO  ,  [ 0 ] , [ ] , T2 , t2 , Period) ; 
p (4) =2 *pi /Period; 
else  %Period==LombPeriod 
p0= [40  .5  1]; 

p=fmins { ' sumsqerrorFP ' , pO , [ 0 ] , [ ] , T2 , t2 , Period) ; 
p (4) =2*pi/Period; 

end 

cosinor .Tpred= (Cos inorTpred(p, t) ) ;  %these  use  original  time  array 
cosinor.t=t; 

cosinor .N=length(t2) ;  %this  is  the  number  of  -isnan  in  T 

cosinor .mesor=p(l) ; 
cosinor .amplitude=p (2) ; 
cosinor .period=2*pi/p (4) ; 

cosinor .phase=-p (3 ) *cosinor. period/2 /pi;  %for  cosine  cosinor 

FirstMidnightAf terTimel=24 ; 

while  FirstMidnightAf terTimel<t ( 1 ) 

FirstMidnightAfterTimel=FirstMidnightAf terTimel+24 ; 

end 

cosinor . gphase=cosinor .phase ; 

while  cosinor . gphase<FirstMidnightAf terTimel 

cosinor . gphase=cosinor . gphase+cosinor .period; 

end 

cosinor . clock_phase=cosinor . gphase ; 
while  cosinor .clock_phase>24 
cosinor . clock_phase=cosinor . clock_phase-24 ; 
end 

cosinor .phase_clockstring=datestr (cosinor .clock_phase/24 , 13) ; 

%correlation  statistics 
Cosinor_Tpred2=CosinorTpred(p, t2) ; 

R=corrcoef (T2 , Cosinor_Tpred2 ) ; 

cosinor.R=R(2) ;  %R  is  a  correlation  matrix.  2nd  coefficient  is  corr  coeff 
%max  and  min  info 

[cosinor .MaxData, iMaxData] =max(T2) ? 

[cosinor .MinData, iMinData] =min(T2) ; 

[cosinor .MaxPred, iMaxPred] =max(Cosinor_Tpred2) ; 

[cosinor .MinPred, iMinPred] =min(Cosinor_Tpred2) ; 
cosinor . TimeOfMaxData=t2 (iMaxData) ; 
while  cosinor .TimeOfMaxData>2 4 
cosinor . TimeOf MaxData=cos inor . TimeOf MaxData-2  4 ; 
end 

cosinor .TimeOf MinData=t2 (iMinData) ; 
while  cosinor .TimeOf MinData>24 
cosinor . TimeOfMinData=cosinor . TimeOfMinData-24  ; 
end 

cosinor  .TimeOf  MaxPred=t2  (iMaxPred)  ; 
while  cosinor .TimeOf MaxPred>2 4 
cosinor . TimeOf MaxPred=cosinor . TimeOf MaxPred-2 4 ; 
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end 

cosinor .TimeOfMinPred=t2 (iMinPred)  ; 
while  cosinor .TimeOfMinPred>24 
cosinor . TimeOfMinPred=cosinor . TimeOfMinPred-24 ; 
end 

return 
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APPENDIX  I 

UTILITY  FUNCTIONS 


% - 

%CircadianDataPath . m 

% - 

function  [Path] =CircadianDataPath 

Path='  j : \ ObservedDa ta \ Circadian\H9 9-11 \ Working  Data\ ' ; 
return 


% - 

%CosinorTpred.m  17  SEP  1999 

% - 

function  Tpred=CosinorTpred (p, t) 

Tpred=p(l) +p(2) *cos (p(4) . *t+p<3) ) ; 
return 


% - 

%DecodeCircadFileName .m 

% - 

function  [Subject ,MPhase, FileNum,MCycle, Julian, Year] =DecodeCircadFileName (CircadFileName) 
f Handle=CircadFileName ( 1 : finds tr ( ' . ' , CircadFileName) -1) ; 

Sub ject=f Handle (2 : 3 ) ; 

MPhase=f Handle (4) ; 

FileNum=str2num(f Handle (5) ) ; 

MCycle=f Handle (6) ; 

Julian=str2num(fHandle(7 : 9 ) )  ; 

Year=str2num (f Handle (11 : 14) )  ; 
return 


% - 

%DecodeJulian.m 

% - 

function  [M, D] =DecodeJulian< Julian, Year) 
cumEOMDay=cumsum (eomday (Year ,1:12)); 
M=min(find{cumEOMDay>=Julian) ) ? 
cumEOMDay= [ 0  cumEOMDay] ; 
D=Julian-cumEOMDay (M) ; 
return 


% - 

%GetFileHandle.m 

% - 

function  [f Handle] =GetFileHandle (fname) 

fHandle=fname (1 : findstr ( ' . ' , fname) -1) ; 
return 


% - 

%InitializeDDE.m 

% - 

function  [channel] =InitializeDDE (topic) 

%initialize  the  file  for  DDE 
channel=ddeinit ( ' excel ' , topic ) ; 

%repeat  the  initialization  step  (bug  in  Matlab  5.3) 
channel=ddeinit ( ' excel ' , topic) ; 
return 


% - 

%PokeTimeBasedColumn.m  11  AUG  1999 


II 


% - - - 

function  PokeTimeBasedColumn (x_heading , x) 
persistent  ColNum 
row=l; 

if  is  empty  (ColNum) 

ColNum=l; 

end 

N=length(x) ; 
if  N>0 

[channel] = Initial izeDDE( 'MATLAB  Results  for  Graphs'); 

Cel l=s treat ( 'r ' ,num2str (row) , 'c' ,num2str (ColNum) ) ; 

ddepoke (channel, Cell, x_heading) ; 
x=ConvertNaNs  (x)  ; 

Cell=strcat ( 'r ' ,num2str (row+1) , 'c' ,num2str (ColNum) , ' :r ' , num2str (N+l) , ' c ' , num2str  (ColNum) ) 

; 

ddepoke (channel , Cell, x) ; 

ColNum=ColNum+l  ; 
ddeterm( channel) ; 

end 

return 

% - 

function  [x] =ConvertNaNs (x) 

MissingValueCode=0; 

IsNotANumber=isnan(x) ; 
index=f ind(IsNotANumber>0) ; 
x( index) =MissingValueCode; 
return 


% - 

%sumsqerror .m  I7  SEP  1999 

% - 

function  sumEsq=sumsqerror (p, T, t) 

sumEsq=sum( (T-CosinorTpred(p, t) ) .*2) ; 
return 


% - 

%sumsqerrorFP.m  I7  SEP  1999 

% - 

function  sumEsq=sumsqerrorFP (p, T, t, Period) 
p (4) =2*pi/Period; 

sumEsq=sum ( (T-CosinorTpred(p, t) ) . A2) ; 
return 


% - 

%TestDDEConnection.m  11  AUG  1999 

% - 

function  TestDDEConnection (Worksheet) 

%Test  that  Excel  file  is  open  and  active 

%InitializeDDE  is  an  external  function  named  InitializeDDE.m 
[channel] =InitializeDDE (Worksheet) ; 
if  channel -=0 

errordlg(strcat  ( 'Excel  worksheet  Worksheet, '' '  must  be  open  and 

active' ) , 'ERROR' ) ; 

error (s treat ( 'Excel  worksheet  ' ' ' , Worksheet, ' ' '  must  be  open  and  active' ) ) ; 

end 

ddeterm ( channel ) ; 
return 


% - 

%WriteSummaryResultsByMethod.m  1  OCT  1999 

% - 

function  WriteSummaryResultsByMethod (FitData) 
persistent  RowNum 
%constants 

Worksheet^ ' MATLAB  Summary  Results'; 

AHeading={ 'mesor' ;  'amplitude' ;  'period' ;  'gphase' ;  ' clock_phase' ; . . . 
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7phase„clockstring7 ;  'R';  7N7 ;  'MaxData7;  7TimeOfMaxData7 ; . . . 

'MinData';  'TimeOfMinData7 ;  7MaxPred7 ;  7TimeOfMaxPred7 ;  'MinPred7 ; 
'TimeOfMinPred7 } ; 

NumericYN= [1  1111011111111X1]; 

%get  row  number  to  start  at 
if  isempty (RowNum) 

RowNum=l; 

end 

startcol=l; 

%open  the  DDE  connection  to  the  excel  worksheet 
[channel] =InitializeDDE (Worksheet) ; 

ndays=length (FitData.R) ; 

ddepoke  (channel,  s treat  { 7r7  ,num2str  (RowNum)  ,  'c'  /num2str  (started)  )  , FitData. Title)  ; 

col=startcol? 

if  ndays>l 

PokeSummaryColumn (channel, RowNum+1, col, ndays, "Day  (rel  to  Event)7,... 

FitData . DayRe la tive2 Event , 1 ) ; col=col+l ; 

else 

col=col+l; 

end 

[NumHeadings, dummy] =size (AHeading) ; 
for  z=l : NumHeadings 

Heading=char (AHeading (z) ) ; 

Data=eval ( streat ( 7  FitData . 7 , Heading) ) ; 

PokeSummaryColumn  ( channel ,  RowNum+1 ,  col ,  ndays ,  Heading ,  Data,NumericYN(z) )  ; 
col=col+l; 

end 

%terminate  DDE  connection 
dde term (channel) ; 

%update  starting  row  number 
RowNum=RowNum+ndays+3 ; 

return 

% - 

function  PokeSummaryColumn (channel , row, col , ndays , x_heading, x, NumericYN) 
ddepoke (channel, streat { 7r 7 , num2str (row) , 7c7 ,num2str (col) ) ,x_heading) ; 
if  NumericYN==0 
for  i=l: ndays 

ddepoke (channel, streat ( 7r 7 , num2str (row+i) , 7c7 ,num2str (col) ) ,x(i, : ) ) ; 

end 

else 

ddepoke (channel, streat (7r 7 ,num2str (row+1) , 7 c 7 ,num2str (col) , . . . 

7 :r7 ,num2str (row+ndays) , 7c7 , num2str (col) ) ,x) ; 

end 


return 
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